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1 The Boltzmann equation 



1.1 The equation 

The Boltzmann equation is the first and the most celebrated nonlinear kinetic equation 
introduced by the great Austrian scientist Ludwig Boltzmann in 1872 £Q. This equation 
describes the dynamics of a moderately rarefied gas, taking into account the two pro- 
cesses: the free flight of the particles, and their collisions. In its original version, the 
Boltzmann equation has been formulated for particles represented by hard spheres. The 
physical condition of rarefaction means that only pair collisions are taken into account, 
a mathematical specification of which is given by the Grad- Boltzmann limit .2]: If iV is 
the number of particles, and a is the diameter of the hard sphere, then the Boltzmann 
equation is expected to hold when N tends to infinity, a tends to zero, Na 3 (the volume 
occupied by the particles) tends to zero, while Na 2 (the total collision cross section) re- 
mains constant. The microscopic state of the gas at time t is described by the one-body 
distribution function P(x,v,t), where x is the position of the center of the particle, and 
v is the velocity of the particle. The distribution function is the probability density of 
finding the particle at time t within the infinitesimal phase space volume centered at the 
phase point (x, v). The collision mechanism of two hard spheres is presented by a relation 
between the velocities of the particles before [v and w] and after [v 1 and w'\ their impact: 

v' = v — n(n, v — w), 
w' = w + n(n, v — w), 

where n is the unit vector along v — v'. Transformation of the velocities conserves the 
total momentum of the pair of colliding particles (v' + w' = v + w), and the total kinetic 
energy (v' 2 + w' 2 = v 2 + w 2 ) The Boltzmann equation reads: 

f + {" d £) = ^/ fl3 / B _(p(*y, f )p ( *v,*) 

—P(x,v,t)P(x,w,t)) I (it; — v,n) | dwdn, (1) 

where integration in w is carried over the whole space R 3 , while integration in n goes over 
a hemisphere B~ = {n G S 2 \ (w—v, n) < 0} . This inequality (w—v, n) < corresponds 
to the particles entering the collision. The nonlinear integral operator in the right hand 
side of (0) is nonlocal in the velocity variable, and local in space. The Boltzmann equation 
for arbitrary hard-core interaction is a generalization of the Boltzmann equation for hard 
spheres under the proviso that the true infinite-range interaction potential between the 
particles is cut-off at some distance. This generalization amounts to a replacement, 

a 2 \(w- v, n) | dn -> B(6, | w - v |) d# de, (2) 

where function B is determined by the interaction potential, and vector n is identified 
with two angles, 6 and e. In particular, for potentials proportional to the n-th inverse 
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power of the distance, the function B reads, 



B(9,\ v-w |) = /9(0) | v-w . (3) 

In the special case n = 5, function B is independent of the magnitude of the relative 
velocity (Maxwell molecules). Maxwell molecules occupy a distinct place in the theory 
of the Boltzmann equation, they provide exact results. Three most important findings 
for the Maxwell molecules are mentioned here: 1. The exact spectrum of the linearized 
Boltzmann collision integral, found by Truesdell and Muncaster 2. Exact transport 
coefficients found by Maxwell even before the Boltzmann equation was formulated, 3. 
Exact solutions to the space-free model version of the nonlinear Boltzmann equation. 
Pivotal results in this domain belong to Galkin [3] who has found the general solution to 
the system of moment equations in a form of a series expansion, to Bobylev, Krook and 
Wu EH [7| who have found an exact solution of a particular elegant closed form, and 
to Bobylev who has demonstrated the complete integrability of this dynamic system [B] , 
the review of relaxation of spatially uniform dilute gases for several types of interaction 
models, of exact solutions and related topics was given in P . 

A broad review of the Boltzmann equation and analysis of analytical solutions to 
kinetic models is presented in the book of Cercignani ^Uj- A modern account of rigorous 
results on the Boltzmann equation is given in the book ^T] . Proof of the existence theorem 
for the Boltzmann equation was given by DiPerna and Lions [T2~] . 

It is customary to write the Boltzmann equation using another normalization of the 
distribution function, f(x, v, t) dx dv, taken in such a way that the function / is compliant 
with the definition of the hydrodynamic fields: the mass density p, the momentum density 
pu, and the energy density e: 

j f(x,v,t)mdv = p{x,t), 

J f(x, v, t)mv dv = pu(x,t), (4) 
f v 2 

/ f(x,v,t)m—dv = e(x,t). 

Here m is the particle's mass. 

The Boltzmann equation for the distribution function / reads, 

(5) 



dt \ dx 

where the nonlinear integral operator in the right hand side is the Boltzmann collision 
integral, 

Q= f [ (f(v')f(w')-f(v)f(w))B(6,v)dwd6d£. (6) 

JRi JB- 
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Finally, we mention the following form of the Boltzmann collision integral (sometimes 
referred to as the scattering or the quasi- chemical representation), 

Q = J W(v,w\v',w')[(f(v')f(w')-f(v)f(w))}dwdw'dv', (7) 

where W is a generalized function which is called the probability density of the elementary 
event, 

W = w(v, w I v', w')S(v + w-v'- w')5{v 2 + w 2 - v' 2 - w' 2 ). (8) 

1.2 The basic properties of the Boltzmann equation 

Generalized function W has the following symmetries: 

W(v',w' | v,w) = W(w',v' | v,w) 

= W(v',w' | w,v) = W(v,w | v',w'). (9) 

The first two identities reflect the symmetry of the collision process with respect to 
labeling the particles, whereas the last identity is the celebrated detailed balance condi- 
tion which is underpinned by the time-reversal symmetry of the microscopic (Newton's) 
equations of motion. The basic properties of the Boltzmann equation are: 

1. Additive invariants of collision operator. 

j Q(f,f){l,v,v 2 }dv = 0, (10) 

for any function /, assuming integrals exist. Equality (|10|) reflects the fact that the number 
of particles, the three components of particle's momentum, and the particle's energy are 
conserved in collisions. Conservation laws (jlUj) imply that the local hydrodynamic fields 
(J3J can change in time only due to redistribution over the space. 

2. Zero point of the integral (Q = 0) satisfy the equation (which is also called the 
detailed balance): For almost all velocities, 

/(«', x, t)f(w', x, t) = f(v, x, t)f(w, X, t). 

3. Boltzmann's local entropy production inequality: 

a(x,t) = -k B J Q(f,f)\nfdv>0, (11) 

for any function /, assuming integrals exist. Dimensional Boltzmann's constant (&b ~ 
6 • 10 _23 Jole/Kelvin) in this expression serves for a recalculation of the energy units into 
the absolute temperature units. Moreover, equality sign takes place if In / is a linear 
combination of the additive invariants of collision. 
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Distribution functions / whose logarithm is a linear combination of additive collision 
invariants, with coefficients dependent on x, are called local Maxwell distribution functions 
/lm, 

/LM = m {-^T ) 6XP { 2k B T ) ■ (12) 

Local Maxwellians are parametrized by values of five hydrodynamic variables, p , u 
and T. This parametrization is consistent with the definitions of the hydrodynamic fields 
®j J fhM{ m , mv , mv 2 /2} dv = (p,pu,e) provided the relation between the energy and 
the kinetic temperature T, holds, s = 2 mk B r - 

4. Boltzmann's H theorem: The function 

S[f] = -k B J flnfdv, (13) 

is called the entropy density 1 . The local H theorem for distribution functions independent 
of space states that the rate of the entropy density increase is equal to the nonnegative 
entropy production, 

f = „>0. (14) 

Thus, if no space dependence is concerned, the Boltzmann equation describes relax- 
ation to the unique global Maxwellian (whose parameters are fixed by initial conditions), 
and the entropy density grows monotonically along the solutions. Mathematical specifi- 
cations of this property has been initialized by Carleman, and many estimations of the 
entropy growth were obtained over the past two decades. In the case of space-dependent 
distribution functions, the local entropy density obeys the entropy balance equation: 

dS(x,t) + /d^ M \ =ff(sM) > 0) (15) 



dt \dx 

where J s is the entropy flux, J s (x, t) = —k B J In f(x, t)vf(x, t) dv. For suitable boundary 
conditions, such as, specularly reflecting or at the infinity, the entropy flux gives no 
contribution to the equation for the total entropy, Stot = J S(x,t)dx and its rate of 
changes is then equal to the nonnegative total entropy production o tot = J a(x, t) dx (the 
global H theorem). For more general boundary conditions which maintain the entropy 
influx the global H theorem needs to be modified. A detailed discussion of this question 
is given by Cercignani ^U]- The local Maxwellian is also specified as the maximizer of the 
Boltzmann entropy function subject to fixed hydrodynamic constraints (jlj). For this 
reason, the local Maxwellian is also termed as the local equilibrium distribution function. 

1 From physical point of view a value of the function / can be treated as dimensional quantity, but 
if one changes the scale and multiplies / by a positive number v then S[f] transforms into vS[f] + 
vhiu J f dv. For a closed system the corresponding transformation of the entropy is an inhomogeneous 
linear transformation with constant coefficients. 
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1.3 Linearized collision integral 



Linearization of the Boltzmann integral around the local equilibrium results in the linear 
integral operator, 



Lh(v) = J W(v,w | v',w')f LM (v)f LM (w) 
h(v') h(w r ) h(v) h(w) 



x 



/lm(«') /lm(w') /lm(w) /lm(w) 



dit/ dv' dw. 



(16) 



Linearized collision integral is symmetric with respect to scalar product defined by the 
second derivative of the entropy functional, 



fLM( v )9(v)Lh(v)dv = / f L ^(v)h(v)Lg(v)dv, 



it is nonpositively definite, 



f L ^(v)h(v)Lh(v) dv<0, 



where equality sign takes place if the function hf^ is a linear combination of collision 
invariants, which characterize the null-space of the operator L. Spectrum of the linearized 
collision integral is well studied in the case of the small angle cut-off. 



2 Phenomenology and Quasi-chemical representation 
of the Boltzmann equation 

Boltzmann's original derivation of his collision integral was based on a phenomenological 
"bookkeeping" of the gain and of the loss of probability density in the collision process. 
This derivation postulates that the rate of gain G equals 

G = J W + (v,w | v',w')f(v')f(w')dv'dw'dw, 

while the rate of loss is 

L = y W~(v, w | v', w')f(v)f(w) dv'dw'dw. 

The form of the gain and of the loss, containing products of one-body distribution 
functions in place of the two-body distribution, constitutes the famous Stosszahlansatz. 
The Boltzmann collision integral follows now as (G — L), subject to the detailed balance 
for the rates of individual collisions, 

W + (v,w | v',w') = W~(v,w | v',w'). 
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This representation for interactions different from hard spheres requires also the cut- 
off of functions /3 (J3J) at small angles. The gain-loss form of the collision integral makes 
it evident that the detailed balance for the rates of individual collisions is sufficient to 
prove the local H theorem. A weaker condition which is also sufficient to establish the H 
theorem was first derived by Stueckelberg ^3] (so-called semi- detailed balance), and later 
generalized to inequalities of concordance 



J dv' J dw'(W + (v, w | v', w') - W~(v, w | v', w')) > 0, 

J dv J dw{W + {v, w | v', w') - W~(v, w | v', w')) < 0. 

The semi-detailed balance follows from these expressions if the inequality signes are 
replaced by equalities. 

The pattern of Boltzmann's phenomenological approach is often used in order to con- 
struct nonlinear kinetic models. In particular, nonlinear equations of chemical kinetics are 
based on this idea: If n chemical species Ai participate in a complex chemical reaction, 

^2a si Ai ^2f3 si Ai, 

i i 

where a S i and (3 s i are nonnegative integers {stoichiometric coefficients) then equations of 
chemical kinetics for the concentrations of species Cj are written 



dt 

s=l 



n dG \ ( " dG 

& ex p i E - ^ ex p l^E dc-/°> 



Functions iff and ip~ are interpreted as constants of the direct and of the inverse 
reactions, while the function G is an analog of the Boltzmann's if-function. 

Modern derivation of the Boltzmann equation, initialized by the seminal work of N.N. 
Bogoliubov Cj\ . seek a replacement condition for the Stosszahlansatz which would be 
more closely related to many-particle dynamics. Different conditions has been formulated 
by D.N. Zubarev ^H]; R.M. Lewis JTj and others. The advantage of these formulations 
is the possibility to systematically find corrections not included in the Stosszahlansatz. 



3 Kinetic models 

Mathematical complications caused by the nonlinearly Boltzmann collision integral are 
traced back to the Stosszahlansatz. Several approaches were developed in order to simplify 
the Boltzmann equation. Such simplifications are termed kinetic models. Various kinetic 
models preserve certain features of the Boltzmann equation, while sacrificing the rest of 
them. The most well known kinetic model which preserves the H theorem is the nonlinear 
Bhatnagar-Gross-Krook model (BGK) 18J. The BGK collision integral reads: 

Qbgk = — (/ — /lm(/))- 
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The time parameter r > is interpreted as a characteristic relaxation time to the local 
Maxwellian. The BGK collision integral is the nonlinear operator: Parameters of the 
local Maxwellian are identified with the values of the corresponding moments of the 
distribution function /. This nonlinearly is of "lower dimension" than in the Boltzmann 
collision integral because /lm(/) is a nonlinear function of only the moments of / whereas 
the Boltzmann collision integral is nonlinear in the distribution function / itself. This 
type of simplification introduced by the BGK approach is closely related to the family 
of so-called mean-field approximations in statistical mechanics. By its construction, the 
BGK collision integral preserves the following three properties of the Boltzmann equation: 
additive invariants of collision, uniqueness of the equilibrium, and the H theorem. A class 
of kinetic models which generalized the BGK model to quasiequilibrium approximations 
of a general form is described as follows: The quasiequilibrium /* for the set of linear 
functionales M(f) is a distribution function f*(M)(x,v) which maximizes the entropy 
under fixed values of functions M. The Quasiequilibrium (QE) models are characterized 
by the collision integral [19J, 



Same as in the case of the BGK collision integral, operator Qqe is nonlinear in the 
moments M only. The QE models preserve the following properties of the Boltzmann 
collision operator: additive invariants, uniqueness of the equilibrium, and the H theorem, 
provided the relaxation time r to the quasiequilibrium is sufficiently small JH] ■ A different 
nonlinear model was proposed by Lebowitz, Frisch and Helfand [20J: 



The collision integral has the form of the self-consistent Fokker-Planck operator, describing 
diffusion (in the velocity space) in the self-consistent potential. Diffusion coefficient D > 
may depend on the distribution function /. Operator Qd preserves the same properties 
of the Boltzmann collision operator as the BGK model. The kinetic BGK model has been 
used for obtaining exact solutions of gasdynamic problems, especially its linearized form 
for stationary problems. Linearized BGK collision model has been extended to model 
more precisely the linearized Boltzmann collision integral [TO] . 

4 Methods of reduced description 

One of the major issues raised by the Boltzmann equation is the problem of the reduced 
description. Equations of hydrodynamics constitute a closed set of equations for the 
hydrodynamic field (local density, local momentum, and local temperature). From the 
standpoint of the Boltzmann equation, these quantities are low-order moments of the one- 
body distribution function, or, in other words, the macroscopic variables. The problem 
of the reduced description consists in giving an answer to the following two questions: 



Qqe(/) = --[f-t(M(f))) + Q(t(M(f)),r(M(f))). 
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1. What are the conditions under which the macroscopic description sets in? 

2. How to derive equations for the macroscopic variables from kinetic equations? 
The classical methods of reduced description for the Boltzmann equation are: the 

Hilbert method, the Chapman- Enskog method, and the Grad moment method. 



4.1 The Hilbert method 

In 1911, David Hilbert introduced the notion of normal solutions, 

fn(v, n(x,t), u(x,t), T(x,t)), 

that is, solutions to the Boltzmann equation which depend on space and time only through 
five hydro dynamic fields [2~T] . 

The normal solutions are found from a singularly perturbed Boltzmann equation, 

D t f = -Q{fJ), (17) 

where e is a small parameter, and 

Physically, parameter e corresponds to the Knudsen number, the ratio between the mean 
free path of the molecules between collisions, and the characteristic scale of variation of the 
hydrodynamic fields. In the Hilbert method, one seeks functions n(x,t), u(x,t), T(x,t), 
such that the normal solution in the form of the Hilbert expansion, 



(0 

c J: 

i=0 



satisfies the f!17|) order by order. Hilbert was able to demonstrate that this is formally 
possible. Substituting (JTHJ) into (fT7|) . and matching various order in e, we obtain the 
sequence of integral equations 

Wh 0) ,/h 0) )=0, (19) 
LfS ] = D t fg\ (20) 
Lfg ) = D t f£ ) -2Q(fg\fU), (21) 

and so on for higher orders. Here L is the linearized collision integral. From (|19jh it follows 
that is the local Maxwellian with parameters not yet determined. The Fredholm 
alternative, as applied to the second equation from (J2Uj) results in 
a) Solvability condition, 



J D t fS ] {l,v,v 2 }dv = 0, 
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which is the set of the compressible Euler equations of the non-viscous hydrodynamics. 
Solution to the Euler equation determine the parameters of the Maxwellian f^. 

b) General solution = fjp 1 + 2 , where is the special solution to the linear 
integral equation (J2Uj) . and /h'* 2 i s Y e t undetermined linear combination of the additive 
invariants of collision. 

c) Solvability condition to the next equation (|2U[) determines coefficients of the function 

(1)2 

/h in terms of solutions to linear hyperbolic differential equations, 

J A(/4 1)1 + /4 1)a ){i.«.« a }d« = o. 

Hilbert was able to demonstrate that this procedure of constructing the normal solution 
can be carried out to arbitrary order n, where the function is determined from the 
solvability condition at the next, (n + l)-th order. In order to summarize, implementation 
of the Hilbert method requires solutions for the functions n(x,t), u(x,t), and T(x,t) 
obtained from a sequence of partial differential equations. 



4.2 The Chapman-Enskog method 

A completely different approach to the reduced description was invented in 1917 by David 
Enskog [22] 5 and independently by Sidney Chapman j2H|. The key innovation was to seek 
an expansion of the time derivatives of the hydrodynamic variables rather than seeking 
the time-space dependencies of these functions, as in the Hilbert method. 

The Chapman-Enskog method starts also with the singularly perturbed Boltzmann 
equation, and with the expansion 



(n) 
CE- 

n=0 



However, the procedure of evaluation of the functions differs from the Hilbert method: 
Q(fSifS) = 0. (22) 

t/S = -wS.Jffl) + ^/S + (».£) & < 23 > 

Operator / dt is defined from the expansion of the right hand side of hydrodynamic 
equation, 

<9 (0) r , r( mv 2 ) f d \ „ (0) , . . 

—{p,pu,e} = - J im,mv,— j ( v,— 1 /^d«. (24) 

From 4221), function is again the local Maxwellian, whereas ()24|) is the Euler equations, 
and d^/dt acts on various functions g(p, pu, e) according to the chain rule, 

0(0) _ Q g 9(0) Qg Q(0) Qg Q(0) 

~dt 9 ~ dp~~dT p + ~dJp~uj~dT {pu) + ~d~e~dt e " 



11 



while the time derivatives of the hydrodynamic fields are expressed using the right 
hand side of (121) . 

The result of the Chapman- Enskog definition of the time derivative is that the 
Fredholm alternative is satisfied by the right hand side of (|23|). Finally, the solution to 
the homogeneous equation is set to be zero by the requirement that the hydrodynamic 
variables as defined by the function /(°) + e/W coincide with the parameters of the local 
Maxwellian /(°); 

J{l,v,v 2 }f$>dv = 0. 
The first correction of the Chapman- Enskog method adds the terms 
—{p,pu,e\ = - <m,mv 

to the time derivatives of the hydrodynamic fields. These terms correspond to the dis- 
sipative hydrodynamics where viscous momentum transfer and heat transfer are in the 
Navier-Stokes and Fourier form. The Chapman-Enskog method was the first true suc- 
cess of the Boltzmann equation since it made it possible to derive macroscopic equation 
without a priori guessing (the generalization of the Boltzmann equation onto mixtures 
predicted existence of the thermo diffusion before it has been found experimentally), and 
to express transport coefficients in terms of microscopic particle's interaction. 

However, higher-order corrections of the Chapman-Enskog method, resulting in hydro- 
dynamic equations with higher derivatives (Burnett hydrodynamic equations) face severe 
difficulties both from the theoretical, as well as from the practical sides. In particular, 
they result in unphysical instabilities of the equilibrium. 



4.3 The Grad moment method 

In 1949, Harold Grad extended the basic assumption of the Hilbert and the Chapman- 
Enskog methods (the space and time dependence of normal solutions is mediated by the 
five hydrodynamic moments) [21]. A physical rationale behind the Grad moment method 
is an assumption of the decomposition of motions: 

(i) . During the time of order r, a set of distinguished moments M' (which include the hy- 
drodynamic moments and a subset of higher-order moments) does not change significantly 
as compared to the rest of the moments M" (the fast dynamics). 

(ii) . Towards the end of the fast evolution, the values of the moments M" become unam- 
biguously determined by the values of the distinguished moments M' . 

(iii) . On the time of order 9 ^> r, dynamics of the distribution function is determined by 
the dynamics of the distinguished moments while the rest of the moments remain to be 
determined by the distinguished moments (the slow evolution period). 
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Implementation of this picture requires an ansatz for the distribution function in order 
to represent the set of states visited in the course of the slow evolution. In Grad's method, 
these representative sets are finite-order truncations of an expansion of the distribution 
functions in terms of Hermite velocity tensors: 

N 

f G (M',v) = f LM (p,u,E,v)[l + Y,CL {a) (M')H {a) (v - «)], (25) 

(«) 

where Ht a \{v — u) are Hermite tensor polynomials, orthogonal with the weight /lm, while 
coefficient a^(M') are known functions of the distinguished moments M', and N is the 
highest order of M' . Other moments are functions of M'\ M" = M"(fa(M')). 

Slow evolution of distinguished moments is found upon substitution of (|2~5j) into 
the Boltzmann equation and finding the moments of the resulting expression (Grad's 
moment equations). Following Grad, this extremely simple approximation can be im- 
proved by extending the list of distinguished moments. The most well known is Grad's 
thirteen-moment approximation where the set of distinguished moments consists of five 
hydrodynamic moments, five components of the traceless stress tensor = J m[(vi — 
u i)( v j ~ u j) ~ $ij(, v ~ ' u ) 2 /3]/di7, and of the three components of the heat flux vector 
Qi = I ( v i - Ui)m(v - u) 2 /2f dv. 

The decomposition of motions hypothesis cannot be evaluated for its validity within 
the framework of Grad's approach. It is not surprising therefore that Grad's methods 
failed to work in situations where it was (unmotivatedly) supposed to, primarily, in 
the phenomena with sharp time-space dependence such as the strong shock wave. On 
the other hand, Grad's method was quite successful for describing transition between 
parabolic and hyperbolic propagation, in particular, the second sound effect in massive 
solids at low temperatures, and, in general, situations slightly deviating from the classi- 
cal Navier-Stokes-Fourier domain. Finally, the Grad method has been important back- 
ground for development of phenomenological nonequilibrium thermodynamics based on 
hyperbolic first-order equation, the so-called EIT (extended irreversible thermodynamics 

4.4 Special approximations 

Special approximations of the solutions to the Boltzmann equation were found for several 
problems, and which perform better than results of "regular" procedures. The most well 
known is the ansatz introduced independently by Mott-Smith and Tamm for the strong 
shock wave problem: The (stationary) distribution function is thought as 

hus(a(x)) = (1 - a(x))f + + a(x)/_, (26) 

where f± are upstream and downstream Maxwell distribution functions, whereas a(x) is 
an undetermined scalar function of the coordinate along the shock tube. 
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Equation for function a(x) has to be found upon substitution of into the Bolltz- 
mann equation, and upon integration with some velocity-dependent function ip(v). Two 
general problems arise with the special approximation thus constructed: Which function 
tp(v) should be taken, and how to find correction to the ansatz like (|26JI ? 

4.5 The method of invariant manifold 

The general problem of reduced description for dissipative system was recognized as the 
problem of finding stable invariant manifolds in the space of distribution functions (2ZII2HI 
123 EH] • The notion of invariant manifold generalizes the normal solution in the Hilbert 
and in the Chapman-Enskog method, and the finite-moment sets of distribution function 
in the Grad method: If Q is a smooth manifold in the space of distribution function, and 
if fa is an element of Q, then Q is invariant with respect to the dynamic system, 



where Tf n Q is the tangent space of the manifold Q at the point fa. Application of 
the invariant manifold idea to dissipative systems is based on iterations, progressively 
improving the initial approximation, and it involves the following steps: construction of 
the thermodynamic projector and iterations for the invariance condition 

4.5.1 Thermodynamic projector 

Given a manifold Q (not obligatory invariant), the macroscopic dynamics on this manifold 
is defined by the macroscopic vector field, which is the result of a projection of vectors 
J (fa) onto the tangent bundle Tfl. The thermodynamic projector Pj n takes advantage 
of dissipativity: 



where DjS |/ n is the differential of the entropy evaluated in fa. 

This condition of thermodynamicity means that each state of the manifold Q is re- 
garded as the result of decomposition of motions occurring near Q: The state fa is the 
maximum entropy state on the set of states fa + kerPj n . Condition of thermodynamicity 
does not define projector completely; rather, it is the condition that should be satisfied by 
any projector used to define the macroscopic vector field, J' n = P? J(fa). For, once the 
condition (|29jl is met, the macroscopic vector field preserves dissipativity of the original 
microscopic vector field J(f): 




(27) 
(28) 



HJ(fa)eT fn Q, for ail / n en, 



kerP; o C kerDfS 



(29) 



D f S \ fn -P; n Wn)) > for all faeQ. 
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The thermodynamic projector is the formalization of the assumption that £7 is the 
manifold of slow motion: If a fast relaxation takes place at least in a neighborhood of Q, 
then the states visited in this process before arriving at /q belong to kerPj n . In general, 
Pi depends in a non-trivial way on /q. 

4.5.2 Iterations for the invariance condition 

The invariance condition for the manifold Q reads, 

Pn(J(fa)) - J(h) = 0, 

here Pq is arbitrary (not obligatory thermodynamic) projector onto the tangent bun- 
dle of fl The invariance condition is considered as an equation which is solved itera- 
tively, starting with initial approximation Qq. On the (n + 1)— th iteration, the correction 
/( n+1 ) = /( n ) + Sf( n+1 ) is found from linear equations, 

DfWfW = P* n J{f {n) ) - J(/ (n) ), 

P*Jf {n+1) = 0, (30) 

here DfJ* is the linear self-adjoint operator with respect to the scalar product by the 
second differential of the entropy DjS |^( n ). 

Together with the above-mentioned principle of thermodynamic projecting, the self- 
adjoint linearization implements the assumption about the decomposition of motions 
around the n'th approximation. The self-adjoint linearization of the Boltzmann colli- 
sion integral Q (J7J) around a distribution function / is given by the formula, 



D f Q SYM 5f = J W{v,w,\v',w 



,v /(«)/(») + /(t0/(«0 



X 



2 

Sf(v') 5f(w') 5f(v) 5f(w) 
[ f(v>) + f{w>) f(v) /(to) 



dw' dv' dw. (31) 



If / = /lm, the self-adjoint operator (pTTjl becomes the linearized collision integral. 
The method of invariant manifold is the iterative process: 

(f (n \p:) - a (n+ V:) - (f (n+1) ,p: + i) 

On the each 1-st step of the iteration, the linear equation (|3T)|) is solved with the projector 
known from the previous iteration. On the each 2-nd step, the projector is updated, 
following the thermodynamic construction. 

The method of invariant manifold can be further simplified if smallness parameters 
are known. 

The proliferation of the procedure in comparison to the Chapman-Enskog method is 
essentially twofold: 
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First, the projector is made dependent on the manifold. This enlarges the set of 
admissible approximations. 

Second, the method is based on iteration rather than a series expansion in a smallness 
parameter. Importance of iteration procedures is well understood in physics, in partic- 
ular, in the renormalization group approach to reducing the description in equilibrium 
statistical mechanics, and in the Kolmogorov- Arnold-Moser theory of finite-dimensional 
Hamiltonian systems. 

4.6 Quasiequilibrium approximations 

Important generalization of the Grad moment method is the concept of the quasiequilib- 
rium approximations already mentioned above (we discuss this approximation in detail in 
a separate section). The quasiequilibrium distribution function for a set of distinguished 
moment M = m(f) maximizes the entropy density S for fixed M. The quasiequilibrium 
manifold Q* (M) is the collection of the quasiequilibrium distribution functions for all ad- 
missible values of M. The quasiequilibrium approximation is the simplest and extremely 
useful (not only in the kinetic theory itself) implementation of the hypothesis about a 
decomposition of motions: If M are considered as slow variables, then states which could 
be visited in the course of rapid motion in the vicinity of Q*(M) belong to the planes 

Fm = {/ | m(f - r(M)) = 0}. 

In that respect, the thermodynamic construction in the method of invariant manifold is a 
generalization of the quasiequilibrium approximation where the given manifold is equipped 
with a quasiequilibrium structure by choosing appropriately the macroscopic variables of 
the slow motion. In contrast to the quasiequilibrium, the macroscopic variables thus 
constructed are not obligatory moments. A textbook example of the quasiequilibrium 
approximation is the generalized Gaussian function for M = {p, pu, P} where Py = 
J ViVjf dv is the pressure tensor. 

The thermodynamic projector P* for a quasiequilibrium approximation was first in- 
troduced by B. Robertson [31] (in a different context of conservative dynamics and for a 
special case of the Gibbs-Shannon entropy). It acts on a function \1/ as follows 

i 

where M = J rriifdv. The quasiequilibrium approximation does not exist if the highest 
order moment is an odd-order polynomial of velocity (therefore, there exists no quasiequi- 
librium for thirteen Grad's moments), and a regularization is then required. Otherwise, 
the Grad moment approximation is the first-order expansion of the quasiequilibrium 
around the local Maxwellian. 
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5 Discrete velocity models 



If the number of microscopic velocities is reduced drastically to only a finite set, the 
resulting discrete velocity models, continuous in time and in space, can still mimic the 
gas-dynamic flows. This idea was introduced in Broadwell's paper in 1963 to mimic the 
strong shock wave [3~2] . 

Further important development of this idea was due to Cabannes and Gatignol in 
the seventies who introduced a systematic class of discrete velocity models [SSI- The 
structure of the collision operators in the discrete velocity models mimics the polynomial 
character of the Boltzmann collision integral. Discrete velocity models are implemented 
numerically by using the natural operator splitting in which each update due to free flight 
is followed by the collision update, the idea which dates back to Grad. One of the most 
important recent results is the proof of convergence of the discrete velocity models with 
pair collisions to the Boltzmann collision integral |34j . 

6 Direct simulation 

Besides the analytical approach, direct numerical simulation of Boltzmann-type nonlinear 
kinetic equations have been developed since mid of 1960's, beginning with the seminal 
works of Bird |3*o^ l3*o] . The basis of the approach is a representation of the Boltzmann 
gas by a set of particles whose dynamics is modeled as a sequence of free propagation and 
collisions. The modeling of collisions uses a random choice of pairs of particles inside the 
cells of the space, and changing the velocities of these pairs in such a way as to comply 
with the conservation laws, and in accordance with the kernel of the Boltzmann collision 
integral. At present, there exists a variety of this scheme known under the common title 
of the Direct Simulation Monte-Carlo method [SSI EHJ • The DSMC, in particular, provides 
data to test various analytical theories. 

7 Lattice Gas and Lattice Boltzmann models 

Since the mid of 1980's, the kinetic-theory based approach to simulation of complex 
macroscopic phenomena such as hydrodynamics has been developed. The main idea of the 
approach is the construction of minimal kinetic system in such a way that their long-time 
and large-scale limit matches the desired macroscopic equations. For this purpose, the 
fully discrete (in time, space, and velocity) nonlinear kinetic equations are considered on 
sufficiently isotropic lattices, where the links represent the discrete velocities of fictitious 
particles. In the earlier version of the lattice methods, the particle-based picture has 
been exploited, subject to the exclusion rule (one or zero particle per lattice link) [the 
lattice gas model |37.j ]• Most of the present versions use the distribution function picture, 
where populations of the links are non-integer [the lattice Boltzmann model [SHI EH EH 
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EJE2|]- Discrete-time dynamics consists of a propagation step where populations are 
transmitted to adjacent links and collision step where populations of the links at each 
node of the lattice are equilibrated by a certain rule. Most of the present versions use the 
BGK-type equilibration, where the local equilibrium is constructed in such a way as to 
match desired macroscopic equations. The lattice Boltzmann method is a useful approach 
for computational fluid dynamics, effectively compliant with parallel architectures. The 
proof of the H theorem for the Lattice gas models is based on the semi-detailed (or 
Stueckelberg's) balance principle. The proof of the H theorem in the framework of the 
lattice Boltzmann method has been only very recently achieved |431 GHJ S3 ESI EZ] • 

8 Minimal Boltzmann models for flows at low Knud- 
sen number 

In this section, we present a new discrete Boltzmann model which has the correct ther- 
modynamics, and which reproduces Navier- Stokes- Fourier equation in the hydrodynamic 
limit. Discrete velocities are taken as zeros of the Hermite polynomials, and the maximum 
entropy Grad moment method, together with the Gauss-Hermite quadrature is used in 
order to derive the discrete H function. Corresponding local equilibria are found analyt- 
ically. The numerical implementation involves a novel discretization scheme consistent 
with the H theorem. Discretization scheme is extended to a derivation of the diffusive 
boundary condition. Some numerical results of the simulation of the model for different 
flow problems are presented. 

Kinetic theory based simulation schemes, in the context of the computational fluid 
dynamics, have attracted considerable attention in the last decade. In particular, the 
lattice Boltzmann method has emerged as an efficient alternative to the traditional Navier- 
Stokes solvers I4()j. In the spirit of the kinetic theory of gases, in this method simple 
pseudo-particle kinetics in introduced in such a way that hydrodynamic equations are 
obtained as its large-scale long-time limit. The method is essentially a particular discrete 
Boltzmann model with a simplified collision mechanism. The emphasis of this model is 
on recovering desired macroscopic dynamics with minimal computational costs. However, 
it turns out that a vital feature of the kinetic theory, the H theorem, was lost in this 
simplification process HH1 HH] • This deficiency in the method had hindered the further 
development of the method as a tool for simulating thermal hydrodynamics ■ 

The aim of this work is to describe a recently proposed method jlZj, which allows 
construction of lattice Boltzmann like methods from the considerations of the classical 
kinetic theory [TTHITH]. The resulting scheme is computationally as efficient as the standard 
lattice Boltzmann method. The scheme is based on preserving continuous (in time) as 
well as the discrete H theorem fHHSj, and hereby it ensures the non- linear stability of 
the resulting numerical algorithm. 
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8.1 Lattice Boltzmann Method 



We shall first describe the standard lattice Boltzmann model. The basic setup for the 
method is as follows: Let /j(x, t) be a population of discrete velocities c$, i = 1, ... ,b, 
at position x and time t. Hydrodynamic fields are the first few moments of populations, 
namely 

b 

c M , c- } = {p, pu a , pDT + pu 2 }, (32) 

i=l 

where p is the mass density of the fluid, pu a is the momentum density, and pDT + pu 2 is 
the energy density. In our notation, a — 1, . . . , D, denotes the spatial directions, where D 
is the spatial dimension. In the case of athermal hydrodynamics, the list of independent 
hydrodynamic fields consists of only the mass and the momentum density. Typically it is 
assumed that during the collision, the populations relax to their equilibrium value with a 
single relaxation time r (Bhatnagar-Gross-Krook (BGK) form of the collision jTUj) . The 
kinetic equation under this assumption reads, 

dtfi + cM^-r-'ifi-f^), (33) 

where the form of the local equilibrium population, / 4 eq remains to be specified. 

The set of discrete velocities is usually chosen to form links of a Bravais lattice (with 
possibly several sub- lattices). This choice of discrete velocity facilitates in an efficient 
discretization using the method of characteristics. The condition that the local hydrody- 
namic variables are invariants of the collision, 

b 

^2f! q {l } Cia, c 2 } = {p, pu a , pDT + pu 2 }, (34) 

i=l 

imposes a set of restrictions on the form of the equilibrium population . (Note that, in the 
athermal case the energy constraint on the local equilibrium is excluded from this list.) 
Chapman-Enskog analysis of (}33|) . with unknown equilibrium satisfying (|34j) . shows that 
in order to recover the desired hydrodynamic equations, in the large-scale long-time limit 
certain higher order moments of the equilibrium, which are involved in the Chapman- 
Enskog analysis, need to be in the same form as in the classical kinetic theory. Further, 
it is assumed that the equilibrium population can be expressed as a polynomial in the 
hydrodynamic fields. These approximation allow for an explicit construction of the lattice 
Boltzmann method with correct hydrodynamics j4*Uj . 

The above mentioned way of constructing kinetic models is not unique for construc- 
tioning the equilibrium or for selection of the set of discrete velocities. It has been long 
known that not every choice of the discrete velocities and the corresponding equilibria 
results in a stable simulation algorithm. In the case of the athermal hydrodynamics, the 
set of the discrete velocities and corresponding polynomials for the equilibrium have been 
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found which lead to a relatively stable realization. However, the problem of stability be- 
comes particularly severe when thermal modes are included. The reason for this failure is 
attributed to the lack of the H theorem for the method. Recently, it was shown how the 
method should be equipped with the H theorem in the isothermal case. The numerical 
stability of the method is considerably enhanced by the incorporation of the H theorem 

8.2 Entropic discrete BGK model 

The key of the construction is the H function which needs to be found for the dynamics 
under consideration |4lH loT] . By definition, the local equilibrium is then the minimizer of 
the corresponding H function [43], under constraints of the local conservation laws (|H4jl . 
Furthermore, if the corresponding local equilibria are known analytically, we can use the 
BGK model. However, an explicit knowledge of the equilibrium is not a fundamental 
requirement for constructing lattice Boltzmann like method. In other cases, when only 
the H function is known analytically, an alternative single relaxation time model can be 
used, which circumvents the problem of finding the local equilibria in a closed form [4*H] . 

The key question is how to find the correct set of discrete velocities and the corre- 
sponding H function. In order to answer this, we remind the reader of the important 
observation on the relation between the discretization of the velocity space and the well 
known Grad's moment method PUHHI]. Namely, if discrete velocities are constructed from 
zeros of the Hermite polynomials, the method of discrete velocity is essentially equivalent 
to Grad's moment method based on the expansion of the distribution function around 
a fixed Maxwellian distribution function. This link ensures that the resulting hydrody- 
namics is the correct one. However, expanding the local equilibria means giving up the 
advantage of having the H theorem j^Sj. This problem is circumvented, if we directly 
evaluate the Boltzmann H function for the discrete case and construct the dynamics using 
this H function. The idea is to use the Gauss-Hermite quadrature in order to construct 
the set of discrete velocities and the H function. This is essentially equivalent to using the 
entropic Grad's method (the maximum entropy approximation) I n the next section 
we describe how the method works. 

8.2.1 Discrete H function 

Boltzmann's H function written in terms of the one-particle distribution function F(x, c) 
is H = J FlnFdc, where c is the continuous velocity. Close to the local equilibrium, 
this integral is naturally approximated using the Gauss-Hermite quadrature. A direct 
evaluation of Boltzmann's H function by a quadrature reads, 




(35) 
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Table 1: Reconstruction of the macroscopic dynamics with the increase in the order of 
the Hermite polynomial. 



Order of 
Polynomial 


Independent 
variables 


Discrete 

velocities in 1-D 


Wei| 


$hts 


Target equation 


2 


P 


±1 


i 

2 


Diffusion 


3 


P, pu a 


0, ±V3T 


2 1 
3' 6 


Athermal Navier-Stokes, 0(u 2 ) 


4 


P, pu a , E 


± a, ±b 


To 
4a 2 ' 


T 

4& 2 


Thermal Navier-Stokes, 0(8 2 ) 
Athermal Navier-Stokes, 0(u 3 ) 



Here W{ is the weight associated with the ith discrete velocity Cj, and the particles mass 
and the Boltzmann constant /cb are set equal to the unity. Discrete discrete function fi(x) 
are related to values of the continuous distribution function at the nodes of the quadrature 
as fi(x) = Wi(2 7rT )( D//2 ) exp(cf/(2 T ))F(x, Cj). The discrete-velocity entropy functions 
(|33|) for various {wi, Cj} is the unique input for all our constructions below. In this setting, 
the set of discrete velocities corresponds to zeroes of the Hermite polynomials in c/y/2T . 
In the next section, it will be shown how the order of the Hermite polynomial, whose 
zeroes constitute the set of discrete velocities, affect the dynamics. 

8.3 Hydrodynamics 

As the order of the Hermite polynomial is increased (this will correspond to increasing 
number of discrete velocities), the discrete H{ WuC .j function becomes a better approxima- 
tion to the continuous Boltzmann H function. Thus, with the increase in the order of the 
Hermite polynomial a better approximation to the hydrodynamics is constructed. This 
behavior is demonstrated in Tabled 

where a = v3 — V6( Tq) 1 ^ 2 , and b = v3 + \/6( Tq) 1 / 2 . In higher dimensions, the 
discrete velocities are products of the discrete velocities in one dimension, and the weights 
are constructed by multiplying the weights associated with the each component direction. 

In order to give a better prospective of the method, we shall discuss a few examples 
of the Table Q in details. 

8.3.1 Athermal hydrodynamics 

If the discrete velocities are formed using roots of the third-order Hermite polynomials (see 
row 1 of table the Navier-Stokes equation is reproduced up to the order 0(u 2 ). The ex- 
plicit expression for the equilibrium distribution, obtained from solving the minimization 
problem, reads: 

/r = pvi x (36) 
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Table 2: Behavior of higher order moments, in comparison to the continuous case. The 
symbol A denotes the difference from the continuous case. 





A P eq 


AO cq 




Athermal case 


0(u A ) 


0(u 3 ) 




Thermal case* 


0(u 8 ) 


0{u6 2 ), 0{u 3 6), and 0{u 5 ) 


0(9 2 ), 0(u 2 9 2 ), and 0{u 4 ) 



(Tq — T)/Tq is the deviation of the temperature from the reference value. 



D 

n 

a=l 



+ (u a /c s )' 



f (2/y^)( Ma /c s ) + v /l + K/c s p X CiaK ^ Cs) 
\ l-w a /(\/3cs) 



where c s = a/TJ^ is speed of sound. Note that the exponent, (c ia / (^/3c s )), in takes 
the values ±1, and only and the resulting expressions for the equilibrium can be sim- 
plified in each dimensions. Equilibria are positive definite for u a < \/3c s . Without 
going into details of derivation, we mention that the equilibrium ()36|) is the product of 
D one-dimensional solutions which were found earlier in the paper [35] via a direct mini- 
mization procedure, and it coincides with the classical equilibrium of the one dimensional 
Broadwell model [5U]. Also, the H function for the two dimensional case was found earlier 
as the solution to functional equation which imposes the requirement of having correct 
stress tensor |51j . The factorization over spatial dimensions is pertinent to the derivation, 
and is quite similar to the familiar property of Maxwell's distribution function. As the 
higher order moments are not enforced, we need to check the behavior of the higher or- 
der moments. Relevant higher-order moments of the equilibrium distribution, needed to 
establish the athermal hydrodynamics in the framework of Chapman-Enskog method are 
the equilibrium pressure tensor, = ^ fi q Ci a Ci/3, the equilibrium third-order moments, 
Q e ap~t = Y^i fi^iaCipc^, and the equilibrium fourth order moment = J2i c iaCi^ 'f* q . 
For the case of the athermal hydrodynamics, to recover Navier-Stokes equation only the 
equilibrium pressure tensor and the equilibrium third-order moments are required to be 
in the correct form. The deviation of these higher order moments, as compared to the 
continuous case is reported in table [21 

We see that the present equilibrium results in the correct pressure tensor only to the 
order 0(u A ). However, this is sufficient for simulations of low Mach number flows. For 
example, for D = 2, even at a high velocity, u a = 0.25\/7o, the error in the pressure is 
less than 2%. 

As the derivation of the equilibrium is based on the entropy function (f55|) which, 
in turn, is produced by the same quadrature as the polynomial equilibria of the lattice 
Boltzmann method jlH]; it is not surprising that expansions of the function to the or- 
der 0(u 2 ) coincides with the polynomial equilibria (for D = 2) of the standard lattice 
Boltzmann method found earlier. 
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8.3.2 Thermal hydrodynamics 

In precisely the same way, the minimal entropic kinetic model for the thermal case requires 
zeroes of fourth-order Hermite polynomials. In order to evaluate Lagrange multipliers in 
the formal solution to the minimization problem, 

f^ = Wl exp (A + B- Cl + Cc 2 t ), 

we first make an important observation that they can be computed exactly for u = and 
any temperature T within the positivity interval, a 2 < T < b 2 : 



(6 2 -a 2 ) a \w b (b 2 -T)J' 

^=^( vSw^ )- Da2c °- (37) 

With this, we find the equilibrium at zero average velocity and arbitrary temperature, 



ow- -i^r (b 2 - T\ (T - a 2 \ \~^ T , 

^n^*n(— ) (— ) • < 38 > 

Factorization over spatial components is clearly seen in this solution. Once the exact 
solution for zero velocity is known, extension to u 7^ is easily found by perturbation. 
The first few terms of the expanded Lagrange multipliers are: 

A = ^- ( r- a 4 2 -r) " 2 + 0( " 4) - 



C - C I ^ 2 -T)-T(b 2 -3T) 

For the actual numerical implementation, the equilibrium distribution function can be 
calculated analytically, up to any order of accuracy required, by this procedure. Accuracy 
of the relevant higher order moments in this case is shown in the Table El Once the error 
in these terms are small, the present model reconstructs the full thermal hydrodynamic 
equations. 

While in the athermal case the closeness of the resulting macroscopic equations to 
the Navier- Stokes equations is controlled solely by the deviations from zero of the average 
velocity, in the thermal regime an additional control is due to variations of the temperature 
away from the reference value. This means that not only the actual velocity should be 
much less than the heat velocity but also the fractional temperature deviation from To 
should be small, |T — T |/T 1. However, by increasing the reference temperature, 
one gets a wider operating window of the present model. Another important remark is 
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about the use of the this thermal model for the athermal Navier-Stokes equation. If the 
temperature is fixed at the reference value T = T , the pressure tensor and the third 
moment Q^jL becomes exact to the order 0(u 5 ), unlike in the second-order accurate 
standard lattice Boltzmann models and the athermal model constructed above. 

8.3.3 Transport coefficients 

When the single relaxation time BGK model (|33|) is used with the present thermal equi- 
librium, the resulting transport coefficients are as follows: For D — 1, the kinematic 
viscosity v is equal to zero, while the thermal conductivity k is k = (3/2)(rp)T. For 
D > 1, we have v = (rp)T, and k = ((D + 2)/2)(rp)T. The equality of Prandtl number, 
CpVK' 1 , to unity is the well-known limitation of the single relaxation time model which 
is cured by using multi-relaxation time models. The density dependence of the transport 
coefficient in the BGK model can be circumvented by renormalizing the relaxation time 
as t' = t p. 

8.4 Thermodynamics 

In our construction of the discrete velocity model, the main focus is on achieving a good 
approximation of the Boltzmann H function. Thus, we can expect that the correct ther- 
modynamics will be preserved (within the accuracy of the discretization), even in the 
discrete case. Indeed, the local equilibrium entropy, S = — kBH{ WuC .}(f eq ), for the ther- 
mal model satisfies the usual expression for the entropy of the ideal monoatomic gas to 
the overall order of approximation of the method, 

S = pk B \n (T D/2 /p) +0(u\6 2 ). (39) 

8.5 Boundary condition 

Once we have developed a systematic way of obtaining the discrete velocity model from the 
continuous case, the same idea can be extended easily to derive the boundary condition. 
The methodology for obtaining the boundary condition is same as that for obtaining 
the discrete H function. We shall start with the boundary condition written in the 
integral form (see JU|) and then use the Gauss-Hermite quadrature to obtain the boundary 
condition for the discrete case. However, the integral appearing in the boundary condition 
is over half space. At this point, we should use a counterpart of Hermite polynomial 
defined on the half space only. However, once we have chosen the quadrature points for 
the interior (for evaluating the H function), we no longer have freedom to choose the 
quadrature points independently for the evaluation of the boundary condition. For that 
reason, we are forced to use the same quadrature in this case too and then we need to 
evaluate the error in each case considered. It turns out that for the case of the diffusive 
wall in the athermal case, this way of evaluation of the half-space integral of the boundary 
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condition is of the accuracy 0(u 2 ), as in the bulk. For the present purpose, a wall OR is 
completely specified at any point (cc € dR) by the knowledge of the inward unit normal 
n, the wall temperature T w and the wall velocity U w . The explicit expression for the 
boundary condition is 

/i=v^ \(f n \\t°*m \/i (^Pw), ((ci-C7 w ) -n> 0), (40) 

Here, £ denotes the molecular velocity in a frame of reference moving with the wall velocity 
and is equals to c — U w . 



8.6 Spatial and time discretization 

In this subsection, we derive a discretization scheme for the discrete Boltzmann equation 
(|33|). which leads to the lattice Boltzmann equation. 

To begin with, we note that the left hand side of (jHSJ) denotes the convection process, 
in which no dissipation is generated (thus no entropy production). On the other hand, the 
right hand side of the equation denotes the generation of the entropy due to the relaxation 
of the population to its equilibrium value (a local event). In order to explore this physical 
picture further, we look at the solution of the ()33|) after time St 



f i (6t)=exp(-5t(c i -V + L))f i (0)+O[[-) , (41) 




where, L denotes the collision operator and in the case of the BGK approximation con- 
sidered here is given as L/j = T _1 (/ 4 cq — /j). After performing an expansion in powers of 
commutator of the collision and the convection, we can write 



f,(St) = <- q >(-StL) [(>xp(-^(c,-- V))/,(0)]+Ol (-) ). (.12) 



V T 



This, expansion has reduced, the problem into two analytically solvable problems of the 
relaxation free convection and of a local relaxation process. The solution reads, 

fi(x,5t) = fi{x-Ci5t,0) (43) 



+ 





" St' 


^1 — exp 






T 



(/r-/«(x-c^t,0)) + O 




The athermal lattice Boltzmann method utilizes this solution in a efficient manner by 
tuning the time step and the grid spacing in such a way that x — CiSt is always a grid 
point. This solves, the convection process in a trivial fashion. In the thermal case, a 
modification in the algorithm is needed to minimize the mismatch of the spatial grid with 
the grid in the velocity space (set of discrete velocities). The aim is chose the St and <5x in 
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such a way that mismatch is reduced to minimum. Presently, we are working to achieve 
such a discretization. 

However, a important modification in (|4~H|) is needed if we want to achieve the zero 
viscosity limit. In the zero viscosity limit (r — > 0) (J4H|) will require 5i — > 0. This problem 
can be avoided if many short collision steps can be lumped together in some ways. This 
can be done by modifying the discrete kinetic equation (J4H|) as 

/i(x, 8t) = /i(x - c^t, 0) + 2f35t (/f - / f (x - Ci 5t, 0)) , (44) 

where the discrete inverse relaxation time is f3 — 1/ (2r + <5t). In the limit of r ^> 5t, (J4U|) 
and the modified ()44j) are identical. However, in the limit of r tending to zero, these two 
equation have two different relaxational behavior. Note, that only short time dynamics is 
changed here to achieve rapid convergence towards hydrodynamics. Here we remind that 
the relation 

St 

~\ ~ 1 

is a Pade approximation rather than a Taylor series expansion. This also explains, why 
in the standard isothermal lattice Boltzmann method the viscosity is v — (1//3 — St)c 2 s /2. 

The discrete equation (jUJ) is successfully used as for various isothermal flow simu- 
lations |4*U] . However, a drawback of this over- relaxation scheme is the loss of H 
theorem. In (jUJ), the positivity of the population is not guaranteed. This may leads to 
numerical instability in many cases (where during the simulation the populations might 
drift far away from the equilibrium) 



cxp 




(45) 



8.6.1 A geometric procedure for restoring the discrete H theorem 

The advantage of the over-relaxation in achieving large time steps can be retained, and 
the problem of the numerical instability can be removed if populations are over-relaxed 
in such a way that H function does not increase in this process jlHl OH] • To do so, we 
modify (jUJ) as 

/i(x, St) = /i(x - Cl St, 0) + aP (if - /i(x - a6t, 0)) , (46) 

where, a is the solution to the equation, 

H(f) = H(f*). (47) 

here / denotes b dimensional vector consisting of all populations and /* = f—a (/ — / eq ). 
First, the distribution is over- relaxed along the path dictated by local collision (with (3 = 
1) to a point of equal entropy. Afterwards, a second collision with the relaxation parameter 
j3 ensures that in the hydrodynamic limit correct viscosity coefficient is recovered. The 
algorithm is illustrated graphically in the Fig. . This way of implementing the H 
theorem in the method ensures the non-linear stability. For the BGK model it can be 
shown that close to local equilibrium a = 25t. The details of the implementation of the 
discrete H theorem is discussed in the paper |47?] . 
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Figure 1: Stabilization procedure. Curves represent entropy levels, surrounding the local 
equilibrium / eq . The solid curve L is the entropy level with the value H(f) = H(f*), 
where / is the initial, and /* is the auxiliary population. The vector A represents the 
collision integral, the sharp angle between A and the vector — Vi? reflects the entropy 
production inequality. The point M is the solution to (|47j) . The result of the collision 
update is represented by the point f(0). The choice of (3 shown corresponds to the 'over- 
relaxation': H(f((3)) > H(M) but H(f(0)) < H(f). The particular case of the BGK 
collision (not shown) would be represented by a vector A B gk, pointing from / towards 
/ eq , in which case M = f q . 

8.7 Numerical Experiments 

We have performed a numerical simulation of the Kramers' problem jl()j . Kramers' prob- 
lem is a limiting case of the plane Couette flow, where one of the plate is moved to infinity, 
while keeping a fixed shear rate. We compare the analytical solution for the slip-velocity 
at the wall calculated for the linearized BGK collision model ^0] with the numerical 
solution in the Fig. |21 

This shows that one important feature of original Boltzmann equation, the Knudsen 
number dependent slip at the wall is retained in the present model. 

In another numerical experiment, the method was tested in the setup of the two- 
dimensional Poiseuille flow. The time evolution of the computed profile as compared to 
the analytical profile obtained from solving the Navier-Stokes equation is demonstrated 
in Fig. El 

8.8 Outlook 

We have constructed the minimal kinetic models for both athermal and thermal hydrody- 
namic simulations. The natural extension of the Gauss- Hermite quadrature discretization 
towards the Boltzmann H function leads to physically sound kinetic models consistent 
with the H theorem. A new approach to the discretization is taken, which ensures non 
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Figure 2: Relative slip observed at the wall in the simulation of the Kramers' problem 
for shear rate a = 0.001, box length L = 32, Voo = a x L = 0.032 (See for details the 
paper 

linear stability of the numerical method. The simulation of the Kramers' problem as a 
function of the Knudsen number demonstrates the usefulness of the present approach for 
the low Knudsen number flows. 




9 Other kinetic equations 

9.1 The Enskog equation for hard spheres 

The Enskog equation for hard spheres is an extension of the Boltzmann equation to 
moderately dense gases. The Enskog equation explicitly takes into account the nonlocality 
of collisions through a two-fold modification of the Boltzmann collision integral: First, the 
one-particle distribution functions are evaluated at the locations of the centers of spheres, 
separated by the nonzero distance at the impact. This makes the collision integral nonlocal 
in space. Second, the equilibrium pair distribution function at the contact of the spheres 
enhances the scattering probability. 

Enskog's collision integral for hard spheres of radius tq is written in the following form 

M- 

Q = [(v - w) ■ n] \x{x, x + r n)f(x,v')f(x + 2r n,w') 

JRt JB- 

— x(a3, x — r n)f(x, v)f(x — 2r n, w)} dw dn, (48) 

where x{ x i y) is the equilibrium pair-correlation function for given temperature and den- 
sity, and integration in w is carried over the whole space R 3 , while integration in n goes 
over a hemisphere B~ = {n e S 2 \ (w — v,n) < 0}. 

The proof of the H theorem for the Enskog equation has posed certain difficulties, and 
has led to a modification of the collision integral p)3| . 
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Figure 3: Development of the velocity profile in the Poiseuille flow. Reduced velocity 
U y (x) = u y /uy max is shown versus the reduced coordinate across the channel x. Solid line: 
Analytic solution. Different lines correspond to different instants of the reduced time 
T = (ut) / (4R 2 ) , increasing from bottom to top. Here, R is the half-width of the channel. 
Symbol: simulation with the present ELBM algorithm. Parameters used are: viscosity 
v = 5.0015 x 10~ 5 (/3 = 0.9997), steady state maximal velocity u Vrn ^ = 1.10217 x 10" 2 . 
Reynolds number Re = 1157. (See for more details in the paper |46j) 



Methods of solution of the Enskog equation are immediate generalizations of those 
developed for the Boltzmann equation, but there is one additional difficulty The Enskog 
collision integral is nonlocal in space. The Chapman-Enskog method, when applied to 
the Enskog equation, is supplemented with a gradient expansion around the homogeneous 
equilibrium state. 



9.2 The Vlasov equation 

The Vlasov equation (or kinetic equation for a self-consistent force) is the nonlinear equa- 
tion for the one-body distribution function, which takes into account a long-range inter- 
action between particles: 

° f+ (°A + ( F °A = , 



dt \ dx J \ dv 

where F = J <&(| x — x' [) \x-x'\ n ( x> ) ^ x> * s ^ e se lf _coris i s teiit force. In this expression 
$(| x — x' \) \x-x'\ ^ s ^ ne m i crosc °pi c force between the two particles, and n(x') is the 
density of particles, defined self-consistently, n(x') = f f(x',v) dv. 

The Vlasov equation is used for a description of collisionless plasmas in which case it 
is completed by the set of Maxwell equation for the electromagnetic field |54j . It is also 
used for a description of the gravitating gas. 

The Vlasov equation is an infinite-dimensional Hamiltonian system jHHj- Many special 
and approximate (wave-like) solutions to the Vlasov equation are known and they describe 
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important physical effects One of the most well known effects is the Landau damping 
~~|: The energy of a volume element dissipates with the rate 



Q 



E 



|2 u{k) d/p 
dv 



k 2 



where fo is the Maxwell distribution function, | E | is the amplitude of the applied 
monochromatic electric field with the frequency u(k) , k is the wave vector. The Landau 
damping is thermodynamically reversible effect, and it is not accompanied with an entropy 
increase. Thermodynamically reversed to the Landau damping is the plasma echo effect. 



9.3 The Fokker-Planck equation 



The Fokker-Planck equation (FPE) is a familiar model in various problems of nonequi- 
librium statistical physics [SHEET We consider the FPE of the form 



dW(x,t) 
dt 



dx 



W 7T U + 7r W 
ox ox 



(49) 



Here W(x, t) is the probability density over the configuration space x, at the time t, while 
U(x) and D(x) are the potential and the positively semi-definite ((y,Dy) > 0) diffusion 
matrix. 

The FPE ()4T)J) is particularly important in studies of polymer solutions jSE EH Ell • Let 
us recall the two properties of the FPE ([4T?j) : (i). Conservation of the total probability: 
/ W(x, t) dx = 1. (ii). Dissipation: The equilibrium distribution, W eq oc exp(— U), is the 
unique stationary solution to the FPE (j4^j) . The entropy, 



S[W) 



W(x,t) In 



W(x,t) 
W eq (x) 



dx, 



(50) 



is a monotonically growing function due to the FPE ()49|). and it arrives at the global 
maximum in the equilibrium. These properties become most elicit when the FPE (I49j) is 
rewritten as follows: 

d t W(x,t) = M w fj}^. , (51) 



sw(x,ty 



where 



Mi 



w 



d_ 
dx 



W(x,t)D(x) 



d_ 
dx 



is a positive semi-definite symmetric operator with kernel 1. The form (|51j) is the dis- 
sipative part of a structure termed GENERIC (the dissipative vector field is a metric 
transform of the entropy gradient) [62^ I63j . 

The entropy does not depend on kinetic constants. It is the same for different details 
of kinetics, and depends only on the equilibrium data. Let us call this property "uni- 
versality". It is known that for the Boltzmann equation there exists only one universal 
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Lyapunov functional. It is the entropy (we do not distinguish functionals which are re- 
lated to reach other by monotonic transformation). But for the FPE there exists a big 
family of universal Lyapunov functionals. Let h(a) be a convex function of one variable 
a > 0, h"{a) > 0, 



S h [W] 



W eq (x)h 



W(x,t) 



dx. 



W eq {x) 

The density of production of the generalized entropy Sh, ®h is nonnegative: 



a h (x) = W eq (x)h" 



W(x,t) 
W eq (x) 



d W(x,t) d W(x,t) 



dx W eq (x 



dx W eq (x) 



> 0. 



(52) 



(53) 



The most important variants for the choice of h: 
h(a) = a In a, Sh is the Boltzmann-Gibbs-Shannon entropy (in the Kullback form [HU 

ESI), 

h(a) = a In a — elna, e > 0, is the maximal family of additive entropies [fifit IH7| EHj 
(these entropies are additive for composition of independent subsystems). 



h(a) 



Si is the family of Tsallis entropies jHSHZOl- These entropies are not additive, 



but become additive after nonlinear monotonous transformation. This property can serve 
as a definition of the Tsallis entropies in the class of generalized entropies 68J. 



10 Equations of chemical kinetics and their reduction 

10.1 Outline of the dissipative reaction kinetics 

We begin with an outline of the reaction kinetics (for details see e. g. the book 71J). Let 
us consider a closed system with n chemical species Ai, . . . , A n , participating in a complex 
reaction. The complex reaction is represented by the following stoichiometric mechanism: 

a fl iAi + . . . + a sn A n ^ ftiAi + . . . + (3 sn A n , (54) 

where the index s = 1, . . . ,r enumerates the reaction steps, and where integers, a si and 
/3 8 i, are stoichiometric coefficients. For each reaction step s, we introduce n-component 
vectors et s and /3 S with components a S i and (3 S i. Notation ~f s stands for the vector with 
integer components = (3 si — a si (the stoichiometric vector). 

For every A{ an extensive variable Ni, "the number of particles of i-th specie", is 
introduced. The concentration of A, is q = Ni/V, where V is the volume. 

Given the stoichiometric mechanism (|54|) . the reaction kinetic equations read: 

r 

N = VJ{c), J(c) = Y,7sW s (c), (55) 

8=1 
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where dot denotes the time derivative, and W s is the reaction rate function of the step s. 
In particular, the mass action law suggests the polynomial form of the reaction rates: 

n n 

W s (c) = W+(c) - W;{c) = fc+(T) Jjcf- - K(T) JJcf- (56) 

8=1 1=1 

where kf(T) and k~{T) are the constants of the direct and of the inverse reactions rates 
of the sth reaction step, T is the temperature. The (generalized) Arrhenius equation gives 
the most popular form of dependence k+(T): 

kf(T) = afT b ° exp(Sf/k B ) exp(-Hf /k B T) } (57) 

where af, bf are constants, Hf are activation ethalpies, Sf are activation entropies. 

The rate constants are not independent. The principle of detailed balance gives the 
following connection between these constants: There exists such a positive vector c eq (T) 
that 

W+ (c eq ) = W~ (c eq ) for all s = 1 , . . . , r. (58) 

The necessary and sufficient conditions for existence of such c eq can be formulate as 
the system of polynomial equalities for {kf}, if the the stoichiometric vectors {7 S } are 
linearly dependent (see, for example, [7T]). 

The reaction kinetic equations (}55|) do not give us a closed system of equations, because 
dynamics of the volume V is not defined yet. Four classical conditions for closure of 
this system are well studied: U, V — const (isolated system, U is the internal energy); 
H, P = const (thermal isolated isobaric system, P is the pressure, H = U + PV is 
the enthalpy), V, T = const (isochoric isothermal conditions); P, T = const (isobaric 
isothermal conditions). For V, T = const we do not need additional equations and data. 
It is possible just to divide equation ()55|) by the constant volume and write 

r 

c = Y,l s W s {c). (59) 

s=l 

For non-isothermal and non-isochoric conditions we do need addition formulae to de- 
rive T and V. For all the four classical conditions the thermodynamic Lyapunov functions 
G for kinetic equations are known: 

U, V — const, Guy = —S/k-Q-, 

V, T — const, Gy t T = F/k B T = U/k B T - S/k B ; 

H, P = const, Gh,p = —S/k B ; 

P,T = const, G P , T = G/k B T = H/k B T - S/k Bl (60) 

where F = U — TS is the free energy (Helmholtz free energy), G = H — TS is the free 
entalphy (Gibbs free energy). All the thermodynamic Lyapunov functions are normalized 
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to the dimensionless scale (if one measures the number of particles in moles, then it 
is necessary to change fee to R). All these functions decrease in time. For classical 
conditions the corresponding thermodynamic Lyapunov functions can be written in the 
form: G, (const, N). The derivatives <9G.(const, N)/dN( are the same functions of c and 
T for all classical conditions: 

, . <9G. (const, N) fi chcmi (c,T) 
^ T) = W = k B T ' (61) 

where /i c hemi(c, T) is the chemical potential of A{. 

Usual G. (const, N) are strictly convex functions of N, and the matrix dfii/dcj is 
positively definite. The dissipation inequality ([6*2~)l holds 

_^ = V( M| J)<0. (62) 

This inequality is the restriction on possible kinetic law and on possible values of kinetic 
constants. 

One of the most important generalization of the mass action law ([5*6*)) is the Marcelin- 
De Donder kinetic function. This generalization [Z21 IZHj is based on ideas of the thermo- 
dynamic theory of affinity 74 . We use the kinetic function suggested in its final form 
in [231 • Within this approach, the functions W s are constructed as follows: For a given 
/a(c,T) ([61)1 ■ and for a given stoichiometric mechanism (|54|) . we define the gain (+) and 
the loss (— ) rates of the sth step, 

w t = ft ex PO, ot a ), W~ = if J exp(/x, (3 S ), (63) 

where <pf > are kinetic factors, (, ) is the standard inner product, the sum of coordinates 
products. 

The Marcelin-De Donder kinetic function reads: W s = — W~, and the right hand 
side of the kinetic equation ([55)1 becomes, 

r 

J = ^7 s {^+exp(/x,a s ) - ipj exp(/x, (3 S )}. (64) 

8=1 

For the Marcelin-De Donder reaction rate ([55)1 . the dissipation inequality ([6*2*)) is particu- 
larly elegant: 

r 

G = ft) - (/i, a.)] {^+e^) - 9 - s e^} < 0. (65) 

s=l 

The kinetic factors should satisfy certain conditions in order to make valid the dissi- 
pation inequality ([6*5)1 . A well known sufficient condition is the detailed balance: 

V>t = <P7i ( 66 ) 
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other sufficient conditions are discussed in detail elsewhere [TH I71| ITS']. 

For ideal systems, function G, is constructed from the thermodynamic data of indi- 
vidual species. It is convenient to start from the isochoric isothermal conditions. The 
Helmholtz free energy for the ideal system is 

F = k B T Ni [In a - 1 + /ioi] + const T ,y , (67) 

i 

where the internal energy is assumed to be a linear function of N in a given interval of c, 
T: 

U = J2 N M T ) = Yl N i( u °i + C viT), 

i i 

where Ui(T) is the internal energy of A\ per particle. It is well known that S = —(dF/dT)v,N 
U = F + TS = F — T(dF/dT) VjN=const , hence, Ui(T) = —k B T 2 dfj, 0i /dT and 

yUoi = $i + u 0i / k B T - (C Vi /k B ) In T, (68) 

where Si = const, Cyi is the A{ heat capacity at constant volume (per particle). 
In concordance with the form of ideal free energy (|67|) the expression for \x is: 

Hi = In a + Si + u i/k B T - (C Vi / k B ) In T. (69) 



For the function fi of the form (|69|). the Marcelin-De Donder equation casts into the 
more familiar mass action law form (J56|) . Taking into account the principle of detailed 
balance (jo^j) we get the ideal rate functions: 

W s (e) = W+{c)-W;{c), 

\y + (c) = (p s (c T)T~^ iasiCv ^ kB e^ iCiai ^ i+Uo ^ kBT ^ Y\ £ 



n 

1=1 

n 



W~(c) = Vs ^T)T-^^- Cv * /kB e^^ +UOl/kBT) JJcf 31 . (70) 



i=l 



where ip s (c,T) is an arbitrary positive function (from thermodynamic point of view). 

Let us discuss further the vector field J(c) in the concentration space (p)9*J). Conser- 
vation laws (balances) impose linear constraints on admissible vectors dc/di: 

(bi, c) = Bi = const, \ b h = 0, % = 1, . . . , I, (71) 

where bi are fixed and linearly independent vectors. Let us denote as B the set of vectors 
which satisfy the conservation laws (fTTj) with given Bi. 

B = {c\{b 1 ,c)=B 1 ,...,{b l ,c) = B l }. 



34 



The natural phase space X of the system (JoTlj) is the intersection of the cone of n- 
dimensional vectors with nonnegative components, with the set B, and dimX = d = n — l. 
In the sequel, we term a vector c G X the state of the system. In addition, we assume 
that each of the conservation laws is supported by each elementary reaction step, that is 

(7.,&i) = 0, (72) 

for each pair of vectors 7 S and bj. 

Reaction kinetic equations describe variations of the states in time. The phase space 
X is positive-invariant of the system (|59)1 : If c(0) G X, then c(t) G X for all the times 
t > 0. 

In the sequel, we assume that the kinetic equation (JoTJj) describes evolution towards 
the unique equilibrium state, c eq , in the interior of the phase space X. Furthermore, we 
assume that there exists a strictly convex function G(c) which decreases monotonically 
in time due to (|SHJ): 

Here VG is the vector of partial derivatives dG/dci, and the convexity assumes that 
the n x n matrices 

H c =\\d 2 G(c)/d Cl dc 3 l (73) 

are positive definite for all c G X. In addition, we assume that the matrices (|73jl are 
invertible if c is taken in the interior of the phase space. 

The function G is the Lyapunov function of the system ([55)1 . and c eq is the point of 
global minimum of the function G in the phase space X. Otherwise stated, the manifold 
of equilibrium states c eq (Bi, . . . , B{) is the solution to the variational problem, 

G -> min for {b h c) = B i7 i = l,...,l. (74) 

For each fixed value of the conserved quantities Bi, the solution is unique. In many 
cases, however, it is convenient to consider the whole equilibrium manifold, keeping the 
conserved quantities as parameters. 

For example, for perfect systems in a constant volume under a constant temperature, 
the Lyapunov function G reads: 

n 

G = 5>PateA$ q )-l]- (75) 

It is important to stress that c eq in (|75|) is an arbitrary equilibrium of the system, 
under arbitrary values of the balances. In order to compute G(c), it is unnecessary to 
calculate the specific equilibrium c eq which corresponds to the initial state c. Let us 
compare the Lyapunov function G (|75jl with the classical formula for the free energy (j67jl . 
This comparison gives a possible choice for c cq : 

lncf = -6 t - u 0t /k B T + (C Vi /k B ) InT. (76) 
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10.2 The problem of reduced description in chemical kinetics 

What does it mean, "to reduce the description of a chemical system"? This means the 
following: 

1. To shorten the list of species. This, in turn, can be achieved in two ways: 

(i) To eliminate inessential components from the list; 

(ii) To lump some of the species into integrated components. 

2. To shorten the list of reactions. This also can be done in several ways: 

(i) To eliminate inessential reactions, those which do not significantly influence the 
reaction process; 

(ii) To assume that some of the reactions "have been already completed", and 
that the equilibrium has been reached along their paths (this leads to dimensional 
reduction because the rate constants of the "completed" reactions are not used 
thereafter, what one needs are equilibrium constants only). 

3. To decompose the motions into fast and slow, into independent (almost-independent) 
and slaved etc. As the result of such a decomposition, the system admits a study 
"in parts". After that, results of this study are combined into a joint picture. 
There are several approaches which fall into this category. The famous method of 
the quasi-steady state (QSS), pioneered by Bodenstein and Semenov, follows the 
Chapman-Enskog method. The partial equilibrium approximations are predeces- 
sors of the Grad method and quasiequilibrium approximations in physical kinetics. 
These two family of methods have different physical backgrounds and mathematical 
forms. 

10.3 Partial equilibrium approximations 

Quasiequilibrium with respect to reactions is constructed as follows: From the list of 
reactions (J54j) . one selects those which are assumed to equilibrate first. Let they be 
indexed with the numbers s\,...,Sk- The quasiequilibrium manifold is defined by the 
system of equations, 

W+ = W~, 1 = 1,...,*. (77) 

This system of equations looks particularly elegant when written in terms of conjugated 
(dual) variables, n = VG: 

( 7si ,M) = 0, z = l,...,fc. (78) 

In terms of conjugated variables, the quasiequilibrium manifold forms a linear subspace. 
This subspace, L x , is the orthogonal completement to the linear envelope of vectors, 
L = lin{7 Sl ,...,7 s J. 
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Quasiequilibrium with respect to species is constructed practically in the same way 
but without selecting the subset of reactions. For a given set of species, A ix , . . . , A ik , 
one assumes that their concentrations evolve fast to the equilibrium, and remain there. 
Formally, this means that in the /c-dimensional subspace of the space of concentrations 
with the coordinates q 1 , . . . , Ci k , one constructs the subspace L which is defined by the 
balance equations, (6j, c) = 0. In terms of the conjugated variables, the quasiequilibrium 
manifold, L x , is defined by equations, 

^ei 1 , (M = 0"i, •••,/%))■ (79) 

The same quasiequilibrium manifold can be also defined with the help of fictitious reac- 
tions: Let g 1; . . . ,g q be a basis in L. Then ()79)1 may be rewritten as follows: 

(^,/*) = 0, i = l,...,q. (80) 

Illustration: Quasiequilibrium with respect to reactions in hydrogen oxidation: Let us 
assume equilibrium with respect to dissociation reactions, H 2 ^ 2H, and, 2 ^ 20, in 
some subdomain of reaction conditions. This gives: 

^1~ C H 2 = ^1 C H> ^2~ c 2 = ^2 c O- 

Quasiequilibrium with respect to species: For the same reaction, let us assume equilibrium 
over H, O, OH, and H 2 02, in a subdomain of reaction conditions. Subspace L is defined 
by balance constraints: 

ch + coh + 2ch 2 o 2 =0, c + c H + 2ch 2 o 2 = 0. 

Subspace L is two-dimensional. Its basis, {gi,g 2 } in the coordinates ch, cq, coh, and 
Ch 2 o 2 reads: 

g x = (1,1,-1,0), g 2 = (2,2,0,-1). 

Corresponding (JHUj) is: 

/^h + f'O — Von, 2/iH + 2/io = /^h 2 o 2 - 

General construction of the quasiequilibrium manifold: In the space of concentration, 
one defines a subspace L which satisfies the balance constraints: 

(bi,L) =0. 

The orthogonal complement of L in the space with coordinates /j, = VG defines then the 
quasiequilibrium manifold fi^. For the actual computations, one requires the inversion 
from n to c. Duality structure /x <-> c is well studied by many authors [7B1 175j. 

Quasiequilibrium projector. It is not sufficient to just derive the manifold, it is also 
required to define a projector which would transform the vector field defined on the space 



37 



of concentrations into a vector field on the manifold. Quasiequilibrium manifold consists 
of points which minimize G on the affine spaces of the form c + L. These affine planes 
are hypothetic planes of fast motions (G is decreasing in the course of the fast motions). 
Therefore, the quasiequilibrium projector maps the whole space of concentrations on 17 ^ 
parallel to L. The vector field is also projected onto the tangent space of f2x parallel to 
L. 

Thus, the quasiequilibrium approximation assumes the decomposition of motions into 
the fast - parallel to L, and the slow - along the quasiequilibrium manifold. In order 
to construct the quasiequilibrium approximation, knowledge of reaction rate constants of 
"fast" reactions is not required (stoichiometric vectors of all these fast reaction are in L, 
7fast £ L, thus, the knowledge of L suffices), one only needs some confidence in that they 
all are sufficiently fast [Sj. The quasiequilibrium manifold itself is constructed based on 
the knowledge of L and of G. Dynamics on the quasiequilibrium manifold is defined as 
the quasiequilibrium projection of the "slow component" of kinetic equations 

10.4 Model equations 

The assumption behind the quasiequilibrium is the hypothesis of the decomposition of 
motions into fast and slow. The quasiequilibrium approximation itself describes slow 
motions. However, sometimes it becomes necessary to restore to the whole system, and 
to take into account the fast motions as well. With this, it is desirable to keep intact one 
of the important advantages of the quasiequilibrium approximation - its independence of 
the rate constants of fast reactions. For this purpose, the detailed fast kinetics is replaced 
by a model equation (single relaxation time approximation) . 

Quasiequilibrium models (QEM) are constructed as follows: For each concentration 
vector c, consider the affine manifold, c + L. Its intersection with the quasiequilibrium 
manifold Ql consists of one point. This point delivers the minimum to G on c + L. Let 
us denote this point as c* L (c). The equation of the quasiequilibrium model reads: 



where r > is the relaxation time of the fast subsystem. Rates of slow reactions are 
computed at the points c* L (c) (the second term in the right hand side of (|8ip. whereas 
the rapid motion is taken into account by a simple relaxational term (the first term in the 
right hand side of (|81|). The most famous model kinetic equation is the BGK equation 
in the theory of the Boltzmann equation ^H] • The general theory of the quasiequilibrium 
models, including proofs of their thermodynamic consistency, was constructed in the paper 



Single relaxation time gradient models (SRTGM) were introduced in the context of 
the lattice Boltzmann method for hydrodynamics pfo*! H5] . These models are aimed at 




(81) 



slow 



ESI- 
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improving the obvious drawback of quasiequilibrium models (|81|): In order to construct 
the QEM, one needs to compute the function, 

c* L (c) = arg min G(x). (82) 

xdc+L, x>0 

This is a convex programming problem. It does not always have a closed-form solution. 

Let g l , . . . ,g k is some orthonormal basis of L. We denote as D(c) the k x k matrix 
with the elements (g { , H c gj), where H c is the matrix of second derivatives of G ()73j) . Let 
C(c) be the inverse of -D(c). The single relaxation time gradient model has the form: 

k 

c = -\ ^(cUg^ VG) + J2^sW s (c). (83) 

i,j=l slow 

The first term drives the system to the minimum of G on c + L, it does not require 
solving the problem (JE2J); an d its spectrum in the quasiequilibrium is the same as in the 
quasiequilibrium model (J%T|) . Note that the slow component is evaluated in the "current" 
state c. 

The first term of equation (j83|) has a simple form 

c=-VadG + yVW s (c), (84) 

slow 

if one calculates the gradient gradG 6 I on the plane of fast motions c + L with the 
entropic scalar product 2 (x,y) = (x,H c y). 

The models (JHTj) and (|%3"j) lift the quasiequilibrium approximation to a kinetic equation 
by approximating the fast dynamics with a single "reaction rate constant" - relaxation 
time r. 



10.5 Quasi-steady state approximation 

The quasi-steady state approximation (QSS) is a tool used in a major number of works. 
Let us split the list of species in two groups: The basic and the intermediate (radicals 
etc). Concentration vectors are denoted accordingly, c s (slow, basic species), and c f (fast, 
intermediate species). The concentration vector c is the direct sum, c = c s © c f . The fast 
subsystem is f!55|) for the component c f at fixed values of cf. If it happens that in this way 
defined fast subsystem relaxes to a stationary state, c f — > c f qss (c s ), then the assumption 
that c f = Cq SS (c) is precisely the QSS assumption. The slow subsystem is the part of 

2 Let us remind that gradG is the Riesz representation of the differential of G in the phase space X: 
G(c + Ac) = G(c) + (gradG(c), Ac) + o(Ac). ft belongs to the tangent space of X and depends on 
the scalar product. From thermodynamic point of view there is only one distinguished scalar product in 
concentration space, the entropic scalar product. Usual definition of gradG as the vector of partial deriva- 
tives corresponds to the standard scalar product (•,•) and to the choice: X is the whole concentration 
space. In equation 184(1 X = c + L and we use the entropic scalar product. 
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the system ()55|) for c s , in the right hand side of which the component c f is replaced with 
Cq SS (c). Thus, J = J s © Jf, where 

c f = J f (c s © c f ), c s = const; c f -> c^c 8 ); (85) 
c s = J s (c s ©< ss (c s )). (86) 

Bifurcations in the system ()85|) under variation of c s as a parameter are confronted to 
kinetic critical phenomena. Studies of more complicated dynamic phenomena in the fast 
subsystem (|8*3|l require various techniques of averaging, stability analysis of the averaged 
quantities etc. 

Various versions of the QSS method are well possible, and are actually used widely, 
for example, the hierarchical QSS method. There, one defines not a single fast subsystem 
but a hierarchy of them, c , . . . , c ffc . Each subsystem c fl is regarded as a slow system 
for all the foregoing subsystems, and it is regarded as a fast subsystem for the following 
members of the hierarchy. Instead of one system of equations (|85j) . a hierarchy of systems 
of lower- dimensional equations is considered, each of these subsystem is easier to study 
analytically. 

Theory of singularly perturbed systems of ordinary differential equations is used to 
provide a mathematical background and refinements of the QSS approximation. In spite 
of a broad literature on this subject, it remains, in general, unclear, what is the smallness 
parameter that separates the intermediate (fast) species from the basic (slow). Reaction 
rate constants cannot be such a parameter (unlike in the case of the quasiequilibrium) . 
Indeed, intermediate species participate in the same reactions, as the basic species (for 
example, H 2 ^ 2H, H + 2 ^ OH + 0). It is therefore incorrect to state that c f evolve 
faster than c s . In the sense of reaction rate constants, c f is not faster. 

For catalytic reactions, it is not difficult to figure out what is the smallness parameter 
that separates the intermediate species from the basic, and which allows to upgrade 
the QSS assumption to a singular perturbation theory rigorously [71]. This smallness 
parameter is the ratio of balances: Intermediate species include a catalyst, and their total 
amount is simply significantly smaler than the amount of all the Cj's. After renormalizing 
to the variables of one order of magnitude, the small parameter appears explicitly. The 
simplest example provides the catalytic reaction A + Z ^ AZ ^ P + Z (here Z is a 
catalyst, A and P are an initial substrate and a product). The kinetic equations are (in 
obvious notations): 

ca = -kf c A c z + k^CAz, 
c z = -kf c A c z + Kcaz + Hcaz - k^czcp, 
caz = kfc A c z - k^CAz - Hcaz + k^czcp, 
c P = kfcAz - k^czcp. (87) 

The constants and the reactions rates are the same for concentrations ca, cp, and for 
cz, caz, and they cannot be a reason for the relative slowness of ca, cp in comparison with 
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cz-i caz- However, there may be another source of slowness. There are two balances for 
this kinetics: ca + cp + caz = Ba, c z + caz = B z . Let us switch to the dimensionless 
variables: 

q A = ca/ Ba, Sp = c P /B A , = c z /B z , q A z = c A z/B Z - 
The kinetic system (j57)l is then rewritten as 



<a = B 2 
<z = B, 



k~ 

-kf^ASz + -fh^z 

DA 



k~ k + 
-kt^ASz + -ft^az + -ft^az ~ kz^zSp 

£>A D A 
By 

SA + SP + —Saz = 1, <iz + Saz = 1; > 0. (88) 

-DA 

For <C i?^ (the total amount of the catalyst is much smaller than the total amount 
of the substrate) the slowness of <;p is evident from these equations (|5S|). 

For usual radicals, the origin of the smallness parameter is quite similar. There are 
much less radicals than the basic species (otherwise, the QSS assumption is inapplicable). 
In the case of radicals, however, the smallness parameter cannot be extracted directly from 
balances Bi (|71|). Instead, one can come up with a thermodynamic estimate: Function 
G decreases in the course of reactions, whereupon we obtain the limiting estimate of 
concentrations of any specie: 

Ci < max Cj, (89) 

G(c)<G(c(0)) 

where c(0) is the initial composition. If the concentration cr of the radical R is small 
both initially and in the equilibrium, then it should remain small also along the path to 
the equilibrium. For example, in the case of ideal G (|75j) under relevant conditions, for 
any t > 0, the following inequality is valid: 

c R [ln(c R (t)/4 q ) - 1] < G(c(0)). (90) 

Inequality (J90"|) provides the simplest (but rather crude) thermodynamic estimate of cr(£) 
in terms of G(c(0)) and uniformly for t > 0. Complete theory of thermodynamic 
estimates of dynamics has been developed in the book |14j . 

One can also do computations without a priori estimations, if one accepts the QSS 
assumption until the values c f stay sufficiently small. It is the simplest way to operate 
with QSS: Just use it until c f are small! 

Let us assume that an a priori estimate has been found, Cj(i) < q max , for each q. 
These estimate may depend on the initial conditions, thermodynamic data etc. With 
these estimates, we are able to renormalize the variables in the kinetic equations (J55j) in 
such a way that the renormalized variables take their values from the unit segment [0, 1]: 
di = Cj/cj max . Then the system (|55|l can be written as follows: 

§ = ^-Mc). (91) 

Lit- L^i may 
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The system of dimensionless parameters, q = q max / maxj q max defines a hierarchy of 
relaxation times, and with its help one can establish various realizations of the QSS 
approximation. The simplest version is the standard QSS assumption: Parameters 6j 
are separated in two groups, the smaller ones, and of the order 1. Accordingly, the 
concentration vector is split into c s © c f . Various hierarchical QSS are possible, with this, 
the problem becomes more tractable analytically 

There exists a variety of ways to introduce the smallness parameter into kinetic equa- 
tions, and one can find applications to each of the realizations. However, the two particular 
realizations remain basic for chemical kinetics: (i) Fast reactions (under a given thermody- 
namic data); (ii) Small concentrations. In the first case, one is led to the quasiequilibrium 
approximation, in the second case - to the classical QSS assumption. Both of these ap- 
proximations allow for hierarchical realizations, those which include not just two but many 
relaxation time scales. Such a multi-scale approach essentially simplifies analytical studies 
of the problem. 



10.6 Thermodynamic criteria for selection of important reac- 
tions 

One of the problems addressed by the sensitivity analysis is the selection of the important 
and discarding the unimportant reactions. In the paper [ZHJ a simple principle was sug- 
gested to compare importance of different reactions according to their contribution to the 
entropy production (or, which is the same, according to their contribution to G). Based 
on this principle, Dimitrov j2H| described domains of parameters in which the reaction 
of hydrogen oxidation, H 2 + O2 + M, proceeds due to different mechanisms. For each 
elementary reaction, he has derived the domain inside which the contribution of this reac- 
tion is essential (nonnegligible). Due to its simplicity, this entropy production principle is 
especially well suited for analysis of complex problems. In particular, recently, a version 
of the entropy production principle was used in the problem of selection of boundary 
conditions for Grad's moment equations [SHI EI]- For ideal systems (J75)l . as well, as for 
the Marcelin-De Donder kinetics the contribution of the sth reaction to G has a 
particularly simple form: 

G s = -W s \ni^-j, G = Y^G S . (92) 

10.7 Opening 

One of the problems to be focused on when studying closed systems is to prepare exten- 
sions of the result for open or driven by flows systems. External flows are usually taken 
into account by additional terms in the kinetic equations (|55|): 

N = VJ{c)+U(c,t). (93) 
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It is important to stress here that the vector field J(c) in equations ()93|) is the same, as 
for the closed system, with thermodynamic restrictions, Lyapunov functions, etc. The 
thermodynamic structures are important for analysis of open systems (|93j) . if the external 
flow II is small in some sense, for example, it is a linear function of c, has small time 
derivatives, etc. There are some general results for such "weakly open" systems, for 
example, the Prigogine minimum entropy production theorem jH2] and the estimations 
of possible of steady states and limit sets for open systems, based on thermodynamic 
functions and stoihiometric equations |14j . 

There are general results for another limit case: for very intensive flow the dynam- 
ics becomes very simple again [21]. Let the flow has a natural structure: U(c,t) = 
Vin{t)c-in(t) ~ v out(t)c(t), where v in and v out are the rates of inflow and outflow, c in {t) is 
the concentration vector for inflow. If v out is sufficiently big, v out (t) > vo for some critical 
value Vq and all t > 0, then for the open system (}9*3"j) the Lyapunov norm exists: for 
any two solutions c 1 (t) and c 2 (t) the function Hc 1 ^) — c 2 (t)|| monotonically decreases in 
time. Such a critical value vq exists for any norm, for example, for usual Euclidian norm 
II -II 2 = (•,•)• 

For arbitrary form of II the system (|9*H|) can loose all signs of being a thermodynamic 
one. Nevertheless, thermodynamic structures often can help in the study of open systems. 

The crucial questions are: What happens with slow/fast motion separation after open- 
ing? Which slow invariant manifolds for the closed system can be deformed to the slow 
invariant manifolds of the open system? Which slow invariant manifold for the closed sys- 
tem can be used as approximate slow invariant manifold for the open system? There exists 
more or less useful technique to seek the answers for specific systems under consideration 

EH ESI 

The way to study an open system as the result of opening a closed system may be 
fruitful. In any case, out of this way we have just a general dynamical system (|93J) and 
no hints what to do with it. 

*** 

Basic introductory textbook on physical kinetics of the Landau and Lifshitz Course 
of Theoretical Physics [S3] contains many further examples and their applications. 

Modern development of kinetics follows the route of specific numerical methods, such 
as direct simulations. An opposite tendency is also clearly observed, and the kinetic theory 
based schemes are increasingly often used for the development of numerical methods and 
models in mechanics of continuous media. 

References 

[1] Boltzmann, L., Lectures on gas theory, University of California Press, 1964. 



43 



[2] Grad, H. Principles of the kinetic theory of gases, in: S. Fliigge, ed., Handbuch der 
Physics, Band 12, Springer, Berlin, 205-294. 

[3] Truesdell, C, Muncaster, R., Fundamentals of Maxwell's Kinetic Theory of a Simple 
Monatomic Gas, Academic Press, NY, 1980. 

[4] Galkin V. S., Kogan M. N., Makashev N. K., Chapman-Enskog generalized method, 
Dokl. ANauk SSSR, 220 (1975), 304-307. 

[5] Bobylev, A. V., Exact solutions to Boltzmann equations, Dokl. Akad. Nauk SSSR 
225 (1975), 1296-1299; One class of invariant solutions to Boltzmann-equation, ibid. 
231 (1976), 571-574. 

[6] Krook, M, Wu, T. T., Formation of Maxwellian tails, Phys Rev Lett, 36 (1976), 
1107-1109. 

[7] Krook, M, Wu, T. T., Exact solutions of Boltzmann-equation, Phys Fluids, 20 (1977), 
1589-1595. 

[8] Bobylev, A. V., Exact-solutions of the nonlinear Boltzmann-equation and the theory 
of relaxation of a Maxwellian gas, Theor. Math. Phys., 60 (1984), 820-841. 

[9] Ernst, M. H., Nonlinear Model-Boltzmann equations and exact solutions, Physics 
Reports, 78 (1981), 1-171. 

[10] Cercignani, C, The Boltzmann equation and its applications, Springer, New York, 
1988. 

[11] Cercignani, C, Illner, R., Pulvirent, M., The mathematical theory of dilute gases, 
Springer, New York, 1994. 

[12] DiPerna, R.J., Lions, P.L. On the Cauchy problem for Boltzmann equation: Global 
existence and weak stability, Ann. Math, 130 (1989), 321-366. 

[13] Stueckelberg E.C.G., Theoreme H et unitarite de S, Helv. Phys. Acta 25, 5 (1952), 
577-580. 

[14] Gorban, A. N., Equilibrium encircling. Equations of chemical kinetics and their ther- 
modynamic analysis, Nauka, Novosibirsk, 1984. 

[15] Bogolyubov, N. N., Dynamic theory problems in statistical physics, Gostekhizdat, 
Moscow, Leningrad, 1946. 

[16] Zubarev, D., Morozov, V., Ropke, G. Statistical mechanics of nonequilibrium pro- 
cesses, V. 1, Basic concepts, kinetic theory, Akademie Verlag, Berlin, 1996, V. 2, 
Relaxation and hydrodynamic processes, Akademie Verlag, Berlin, 1997. 



44 



[17] Lewis, R. M., A unifying principle in statistical mechanics, J. Math. Phys., 8 (1967), 
1448-1460. 

[18] Bhatnagar, P. L., Gross, E. P., Krook, M., A model for collision processes in gases. 
I. Small amplitude processes in charged and neutral one-component systems, Phys. 
Rev., 94, 3 (1954), 511-525. 

[19] Gorban, A. N., Karlin, I. V., General approach to constructing models of the Boltz- 
mann equation, Physica A, 206 (1994), 401-420. 

[20] Lebowitz, J., Frisch, H., Helfand, E., Non-equilibrium distribution functions in a 
fluid, Physics of Fluids, 3 (1960), 325. 

[21] Hilbert, D. Begriindung der kinetischen Gastheorie, Mathematische Annalen, 72 
(1912), 562-577. 

[22] Enskog, D., Kinetische theorie der Vorange in massig verdunnten Gasen. I Allge- 
meiner Teil, Almqvist and Wiksell, Uppsala, 1917. 

[23] Chapman, S., Cowling, T., Mathematical theory of non-uniform gases, Third edition, 
Cambridge University Press, Cambridge, 1970. 

[24] Grad, H. On the kinetic theory of rarefied gases, Comm. Pure and Appl. Math. 2 4, 
(1949), 331-407. 

[25] Jou, D., Casas- Vazquez, J., Lebon, G., Extended irreversible thermodynamics, 
Springer, Berlin, 1993. 

[26] Muller, I., Ruggeri, T., Extended Thermodynamics, Springer, NY, 1993. 

[27] Gorban, A. N., Karlin, I. V., Zinovyev, A. Yu., Constructive methods of invariant 
manifolds for kinetic problems, Phys. Reports, 396, 4-6 (2004), 197-403. Preprint 
online: http: / / arxiv.or g/abs/cond-mat/0311017. 



[28] Gorban, A. N., Karlin, I. V., The constructing of invariant manifolds for the Boltz- 
mann equation, Adv. Model, and Analysis C, 33(3) (1992), 39-54. 

[29] Gorban, A. N., Karlin, I. V., Thermodynamic parameterization, Physica A, 190 
(1992), 393-404. 

[30] Gorban, A. N., Karlin, I. V., Method of invariant manifolds and regularization of 
acoustic spectra, Transport Theory and Stat. Phys., 23 (1994), 559-632. 

[31] Robertson, B., Equations of motion in nonequilibrium statistical mechanics, Phys. 
Rev., 144 (1966), 151-161. 



45 



[32] Broadwell, J. E., Study of shear flow by the discrete velocity method, J. Fluid Mech. 
19 (1964), 401-414. 

[33] Gatignol, R. Theorie cinetique des gaz a repartition discrete de vitesses. Lecture 
notes in physics, V.36, Springer, Berlin, etc, 1975. 

[34] Palczewski, A., Schneider, J., Bobylev, A.V., A consistency result for a discrete- 
velocity model of the Boltzmann equation, SIAM Journal on Numerical Analysis, 
34, 5 (1997), 1865-1883. 

[35] Bird G.A. Molecular gas dynamics and the direct simulation of gas flows, Clarendon 
Press, Oxford, 1994. 

[36] Oran, E. S. and Oh, C. K. and Cybyk, B. Z., Direct simulation Monte Carlo: recent 
advances and applications, Annu Rev. Fluid Mech., 30 (1998), 403-441. 

[37] Frisch, U., Hasslacher, B. and Pomeau, Y., Lattice-gas automata for the Navier- 
Stokes equation, Phys. Rev. Lett., 56 (1986), 1505-1509. 

[38] Mcnamara, Gr., Zanetti, G., Use of the Boltzmann-equation to simulate lattice-gas 
automata, Phys. Rev. Lett., 61 (1988), 2332-2335. 

[39] Chen S. and Doolen G.D. Lattice Boltzmann method for fluid flows, Annu. Rev. 
Fluid. Mech. 30 (1998), 329-364. 

[40] Succi S., The lattice Boltzmann equation for fluid dynamics and beyond, Clarendon 
Press, Oxford, 2001. 

[41] Higuera, F., Succi, S., Benzi, R., Lattice gas - dynamics with enhanced collisions, 
Europhys. Lett., 9 (1989), 345-349. 

[42] Benzi, R., Succi, S., Vergassola, M., The lattice Boltzmann-equation - theory and 
applications Physics Reports, 222, 3 (1992), 145-197. 

[43] Succi S., Karlin I. V., Chen H., Role of the H theorem in lattice Boltzmann hydro- 
dynamic simulations, Rev. Mod. Phys., 74 (2002), 1203-1220. 

[44] Karlin, I. V., Gorban, A. N., Succi, S., Boffi, V., Maximum entropy principle for 
lattice kinetic equations, Phys. Rev. Lett., 81 (1998), 6-9. 

[45] Ansumali, S., Karlin, I. V., Stabilization of the Lattice Boltzmann method by the H 
theorem: A numerical test, Phys. Rev. E, 62 (6), (2000), 7999-8003. 

[46] Ansumali, S., Karlin, I. V., Entropy function approach to the lattice Boltzmann 
method, J. Stat. Phys., 107 (1/2) (2002), 291-308. 



46 



[47] Ansumali, S., Karlin, I. V., Ottinger, H. C, Minimal entropic kinetic models for 
hydrodynamics, Europhys. Lett., 63 (2003), 798-804. 

[48] Ansumali S., Karlin, I. V. Single relaxation time model for entropic Lattice Boltz- 
mann methods, Phys. Rev. E, 65 (2002), 056312. 

[49] Shan, X., He, X., Discretization of the velocity space in the solution of the Boltzmann 
equation, Phys. Rev. Lett., 80 (1998), 65-67. 

[50] Broadwell, J. E., Shock structure in a simple discerte velocity gas, Phys. Fluids, 7 
(1964), 1243-1247. 

[51] Karlin, I. V., Ferrante, A., Ottinger, H. C. Perfect entropy functions of the Lattice 
Boltzmann method, Europhys. Lett., 47 (1999), 182-188. 

[52] Ansumali, S., Karlin, I. V., Kinetic Boundary condition for the lattice Boltzmann 
method, Phys. Rev. E, 66 (2002), 026311. 

[53] Van Beijeren, H, Ernst, M.H., Modified Enskog equation, Physica A, 68, 3 (1973), 
437-456. 

[54] Lifshitz E.M. and Pitaevskii L.P., Physical kinetics (Landau L.D. and Lifshitz E.M. 
Course of Theoretical Physics, V.10), Pergamon Press, Oxford, 1968. 

[55] Marsden, J.E., Weinstein, A., The Hamiltonian structure of the Maxwell- Vlasov 
equations, Physica D, 4 (1982), 394-406. 

[56] Braun W, Hepp K, Vlasov dynamics and its fluctuations in 1-N limit of interacting 
classical particles, Comm. Math. Phys., 56, 2 (1977), 101-113. 

[57] Van Kampen, N. C, Stochastic processes in physics and chemistry, North-Holland, 
Amsterdam 1981. 

[58] Risken, H., The Fokker-Planck equation, Springer, Berlin, 1984. 

[59] Bird, R. B., Curtiss, C. F., Armstrong, R. C, Hassager, O., Dynamics of Polymer 
Liquids, 2nd edn., Wiley, New York, 1987. 

[60] Doi, M., Edwards, S. F., The theory of polymer dynamics, Clarendon Press, Oxford, 
1986. 

[61] Ottinger, H. C, Stochastic processes in polymeric fluids, Springer, Berlin, 1996. 

[62] Grmela, M., Ottinger, H. C, Dynamics and thermodynamics of complex fluids. I. 
Development of a general formalism, Phys. Rev. E 56 (1997), 6620-6632. 



47 



[63] Ottinger, H. C, Grmela, M., Dynamics and thermodynamics of complex fluids. II. 
Illustrations of a general formalism, Phys. Rev. E, 56 (1997), 6633-6655. 

[64] Kullback, S., Information theory and statistics, Wiley, New York, 1959. 

[65] Plastino, A.R., Miller, H.G., Plastino, A., Minimum Kullback entropy approach to 
the Fokker-Planck equation, Physical Review E 56 (1997). 3927-3934. 

[66] Gorban, A. N., Karlin, I. V., Family of additive entropy functions out 
of thermodynamic limit, Phys. Rev. E, 67 (2003), 016104. Preprint online: 
|http:/ /arxiv.org/abs/cond-ma t / 0205 5TT) 

[67] Gorban, A. N., Karlin, I. V., Ottinger H. C, The additive generalization 
of the Boltzmann entropy, Phys. Rev. E, 67, 067104 (2003). Preprint online: 
http: / / arxiv.org/abs/cond-mat / 0209319 , 

[68] Gorban, P., Monotonically equivalent entropies and solution of ad- 
ditivity equation, Physica A, 328 (2003), 380-390. Preprint online: 
http:// arxiv.ore; /pdf/cond-mat/0304131, 

[69] Tsallis, C, Possible generalization of Boltzmann-Gibbs statistics. J. Stat. Phys., 52 
(1988), 479-487. 

[70] Abe, S., Okamoto, Y. (Eds.), Nonextensive statistical mechanics and its applications, 
Springer, Heidelberg, 2001. 

[71] Yablonskii, G. S., Bykov, V. I., Gorban, A. N., Elokhin, V. I., Kinetic models of 
catalytic reactions. Comprehensive Chemical Kinetics, Vol. 32, Compton R. G. ed., 
Elsevier, Amsterdam (1991). 

[72] Feinberg, M., Chemical kinetics of a sertain class, Arch. Rat. Mech. Anal., 46, 1 
(1972), 1-41. 

[73] Bykov, V. I., Gorban, A. N., Yablonskii, G. S., Description of nonisothermal reactions 
in terms of Marcelin - de Donder kinetics and its generalizations, React. Kinet. Catal. 
Lett., 20, 3-4 (1982), 261-265. 

[74] De Donder, T., Van Rysselberghe, P., Thermodynamic theory of affinity A book of 
principles. Stanford: University Press, 1936. 

[75] Dukek, G., Karlin, I. V., Nonnenmacher, T. F., Dissipative brackets as a tool for 
kinetic modeling, Physica A, 239(4) (1997), 493-508. 

[76] Orlov, N. N., Rozonoer, L. I., The macrodynamics of open systems and the variational 
principle of the local potential, J. Franklin Inst., 318 (1984), 283-314 and 315-347. 



48 



[77] Volpert, A. I., Hudjaev, S. I. Analysis in classes of discontinuous functions and the 
equations of mathematical physics. Dordrecht: Nijhoff, 1985. 

[78] Bykov, V. I., Yablonskii, G. S., Akramov, T. A., The rate of the free energy decrease 
in the course of the complex chemical reaction. Dokl. Akad. Nauk USSR, 234, 3 
(1977) 621-634. 

[79] Dimitrov, V.I., Simple kinetics, Nauka, Novosibirsk, 1982. 

[80] Struchtrup, H., Weiss, W. Maximum of the local entropy production becomes mini- 
mal in stationary processes, Phys. Rev. Lett., 80 (1998), 5048-5051. 

[81] Grmela, M., Karlin, I. V., Zmievski, V. B., Boundary layer minimum entropy prin- 
ciples: A case study, Phys. Rev. E, 66 (2002), 011201. 

[82] Prigogine, I., Thermodynamics of irreversible processes, Interscience, New York, 
1961. 

[83] Zmievskii, V. B., Karlin, I. V., Deville, M., The universal limit in dynamics of dilute 
polymeric solutions, Physica A, 275(1-2) (2000), 152-177. 



49 



